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Abstract  We  report  the  main  results  produced  within  the  effort  sponsored  by  the  U.S.  Air  Force 
Research  Laboratory  with  Contract  No.  FA8718-09-C-0061.  We  focus  on  the  most  relevant 
aspects  of  our  findings,  which  we  have  fully  addressed  during  our  effort:  guided  propagation  and 
leaky-wave  radiation  along  linear  arrays  of  nanoparticles,  also  considering  and  modeling  the 
realistic  presence  of  technological  disorder,  comparison  of  the  guidance  properties  along  linear 
and  planar  arrays  of  nanoparticles  and  nanovoids  in  different  realistic  geometries,  guided  and 
leaky  modes  along  parallel  arrays  of  nanoparticles,  propagation  along  periodic  arrays  of 
nanoparticles  and  their  rigorous  homogenization  as  metamaterials  for  a  variety  of  applications  of 
interest  to  the  U.S.  Air  Force. 
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1.  Executive  Summary 

In  the  present  effort  we  have  explored  a  variety  of  theoretical  and  numerical  problems  associated 
with  wave  propagation  and  radiation  along  periodic  arrays  of  nanoparticles.  Our  extensive  efforts 
have  tackled  several  relevant  problems  in  the  following  areas:  (a)  leaky-  and  guided  modes  along 
isolated  and  coupled  arrays  of  metamaterial  resonators  and  plasmonic  nanoparticles;  (b) 
quantification  of  the  effect  of  degree  of  disorder  in  the  propagation  along  linear  arrays;  (c) 
propagation  along  planar  arrays  of  nanoparticles  and  voids  in  plasmonic  substrates,  and 
comparisons  among  different  solutions  for  nanoscale  light  propagation  and  radiation;  (d) 
rigorous  modeling  of  the  propagation  along  3-D  arrays  of  nanoparticles  and  their  proper 
homogenization  models.  These  results  may  have  a  variety  of  important  applications  of  interest  to 
the  U.S.  Air  Force,  in  particular  in  the  realization  and  design  of  sub-wavelength  compact  optical 
waveguides  with  optimized  field  confinement  and  propagation  length,  nanoantenna  and  leaky- 
wave  radiation  applications  with  optimized  directivity,  control  of  propagation  and  scattering 
from  nanoscale  structures  and  design  of  metamaterials  with  unconventional  wave  interaction.  We 
have  analyzed  several  geometries  and  designs  with  a  variety  of  analytical  and  numerical 
methods,  including  full-wave  numerical  techniques,  modeling  in  details  the  properties  of  the 
arrays  in  terms  of  frequency  and  spatial  dispersion  and  of  radiation  and  Ohmic  losses.  During 
this  effort,  we  have  developed  novel  analytical,  theoretical  and  numerical  tools  for  a  thorough 
understanding  of  the  propagation  and  radiation  properties  of  arrays  of  nanoresonators  in  various 
frequency  regimes  of  operation,  even  in  the  presence  of  realistic  disorder.  Our  efforts  have  been 
successful  not  only  from  the  research  perspective,  but  we  have  also  involved  a  few  students  and 
one  postdoctoral  researcher  within  this  exciting  project.  In  the  following,  we  review  some  of  the 
most  relevant  findings  obtained  in  this  work,  as  outlined  in  recent  publications  referred  at  the  end 
of  this  report.  Each  section  will  contain  one  of  the  topics  developed  in  the  framework  of  this 
project  and  the  main  results  obtained  in  this  framework.  The  interested  reader  is  referred  to  the 
series  of  publications  that  have  been  produced  during  this  effort,  which  are  reported  in  Section  9 
of  this  report. 


2.  Methods,  Assumptions  and  Procedures 

In  this  effort  we  have  used  a  variety  of  analytical  techniques,  including  periodic  Green’s 
functions  and  analytical  continuation  in  the  complex  domain,  as  well  as  numerical  techniques, 
including  finite-integration  technique  and  finite-time-domain  simulations.  We  detail  more  of  our 
numerical  approaches  and  efforts  in  the  following  sections,  specifying  their  relevance  in  each 
part  of  our  effort.  We  assume  in  the  following  sections  an  e~101t  time  dependence  for  the 
electromagnetic  fields,  unless  otherwise  noted. 


3.  Propagation  Along  Linear  Arrays  of  Nanoparticles:  Effect  of  Disorder 

a.  Summary 

Metallic  nanoparticles  are  known  to  possess  several  exciting  properties  when  interacting  with 
light.  Linear  arrays  of  plasmonic  nanoparticles,  in  particular,  can  support  optical  beam 
propagation  with  ultra-confined  lateral  cross-section  (possibly  much  smaller  than  the  free-space 
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wavelength  of  interest),  surpassing  by  far  the  transverse  size  of  regular  optical  fibers.  Although 
in  theory  these  waveguides  may  support  relatively  long  propagation,  in  practical  scenarios  short 
propagation  lengths  have  been  measured.  Absorption  in  plasmonic  materials  indeed  plays  a 
relevant  role  in  this  regard,  but  fabrication  disorder  has  been  also  suspected  to  severely  affect  the 
overall  propagation  length. 

During  this  effort,  we  have  been  able  to  quantify  in  a  simple  analytical  expression  the  role  of 
small  imperfections  in  linear  arrays  of  plasmonic  nanoparticles,  showing  that  the  dominant  effect 
of  random  Gaussian  disorder  on  the  propagation  properties  of  these  arrays  is  equivalent  to  an 
additional  loss  tenn  in  the  metal.  The  total  measured  loss  may  therefore  be  associated  with 
material  absorption  and  to  an  additional  factor  that  is  directly  proportional  to  the  disorder 
variance  and  function  of  the  phase  velocity  of  the  guided  mode.  Our  results  may  not  only  affect 
the  understanding  of  plasmonic  nanoscale  waveguides  and  quantify  the  role  of  technological 
disorder  in  optical  communications,  but  they  may  also  be  extended  to  other  periodic  optical 
devices  of  interest,  like  photonic  crystals  and  metamaterials,  whose  perfonnance  are  often 
limited  by  inherent  fabrication  imperfections.  Analogous  concepts  may  be  able  to  quantify  their 
effect  in  various  practical  designs  of  interest. 

b.  Introduction 

Exploring  the  science  and  applications  of  photonic  crystals  and  metamaterials  formed  by 
periodic  collections  of  inclusions  has  rapidly  grown  in  interest  and  importance  in  recent  years 
(see  e.g.,  [1]).  With  a  few  exceptions  (see,  e.g.,  [2]-[7]),  however,  in  majority  of  the  theoretical 
contributions  on  this  subject  the  role  of  disorders  and  imperfections  due  to  inherent  fabrication 
limitations  are  neglected,  while,  at  the  present  state  of  the  art,  experimental  efforts  in 
metamaterial  realization  and  measurement  inevitably  involve  such  problems.  The  sub¬ 
wavelength  scale  in  which  the  constituent  inclusions  of  a  metamaterial  are  typically  embedded  in 
a  host  medium  usually  allows  a  homogenization  procedure,  suggesting  that  the  intrinsic 
deviation  from  an  ideal  periodic  arrangement  may  not  play  a  fundamental  role  in  the 
metamaterial  perfonnance.  This  disorder,  however,  should  be  taken  into  account  for  a  complete 
and  rigorous  analysis  of  the  metamaterial  response,  which  may  describe  more  closely  the 
experimental  conditions.  A  clear  example  of  such  limitations  is  seen  when  dealing  with  periodic 
lattices  of  inclusions  and  metamaterials:  it  is  well  known  that  an  infinite  periodic  arrangement  of 
lossless  inclusions  cancels  out  the  scattering  losses  from  each  individual  particle,  ensuring  a 
completely  lossless  bulk  metamaterial  [8],  but  the  measured  losses  in  metamaterials,  manifested 
in  the  imaginary  parts  of  extracted  permittivity  and  permeability,  are  often  higher  than  expected 
from  purely  theoretical  considerations.  The  role  of  fabrication  imperfections  and  random 
disorders  enters  directly  into  play  in  this  scenario,  generating  incoherent  superposition  of 
scattered  waves  from  the  different  inclusions,  which  results  in  energy  decay  for  a  wave  traveling 
through  the  metamaterial  sample.  This  energy  “loss”  has  been  often  interpreted  in  the 
measurements  as  energy  absorbed  in  the  bulk  sample,  even  though  it  is  actually  related  mostly  to 
incoherent  scattering,  rather  than  material  absorption. 

A  different,  but  related  problem  that  we  have  recently  analyzed  is  represented  by  the  analytical 
solutions  for  the  eigenmodes  supported  by  1-D  or  3-D  periodic  lattices  of  closely-packed 
plasmonic  nanoparticles  [9]-[  1 0],  which  may  support  highly-confined  guided  modes  and 
backward-wave  propagation,  with  various  potential  applications  for  negative-index  materials  and 
laterally-confined  waveguides.  Also  in  this  case,  the  chain  periodicity  allows  total  cancellation  of 
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the  scattering  losses,  supporting  the  propagation  of  low-loss  eigenmodes  (with  loss  only  due  to 
the  material  absorption  in  the  particles).  However,  disorder  is  inherently  associated  with  the 
realization  of  these  chains:  techniques  for  manufacturing  large  arrays  of  nanoparticles  are 
currently  available,  but  while  they  can  ensure  a  good  control  on  the  average  dimension  of  the 
particles,  it  may  be  more  challenging  to  ensure  a  perfect  alignment  and  positioning  of  each  of 
them  on  a  large  scale.  Since  the  assumption  of  ideally  periodic  arrangement  of  the  particles  has 
provided  us  with  closed-form  solutions  for  the  propagation  of  such  modes  [9]-[  1 0],  which  may 
also  easily  take  into  account  the  presence  of  Ohmic  absorption,  we  may  ask  whether  an  extension 
of  these  analytical  results  may  also  model  the  intrinsic  disorder  expected  in  the  fabrication 
process.  More  broadly,  this  may  be  applied  to  the  general  analysis  of  metamaterials,  in  order  to 
quantitatively  model  how  the  deviation  from  ideal  periodic  arrangements  of  the  inclusions  may 
affect  their  potential  application  in  realistic  setups  and  devices. 

In  this  section,  we  explore  these  concepts  and  in  particular  we  study  how  the  presence  of 
disorder  may  be  modeled  in  simple  geometries  like  periodically  arranged  nanoparticles.  In 
particular,  we  show  how  such  disorder  effects,  when  small  enough,  may  be  equivalently 
described  as  the  effect  of  an  additional  “loss”  in  the  materials  forming  the  particles,  directly 
related  to  the  statistical  variance  of  the  disorder  in  the  system,  and  may  characterize  both 
quantitatively  and  qualitatively  the  effects  of  imperfections  in  the  realization  of  these 
metamaterials.  These  results  may  be  further  generalized  to  generic  metamaterial  lattices, 
providing  relevant  physical  insights  into  the  analysis  of  metamaterials  for  a  practical  realization. 

It  should  be  pointed  out  that  other  groups  have  recently  analyzed  the  problem  of  disorders  in 
metamaterials,  both  from  theoretical  and  experimental  points  of  view  (see  e.g.,  [2]-[5]),  showing 
how  these  effects,  although  of  second-order,  may  be  relevant  in  some  applications  and  how  they 
are  sometimes  underestimated.  An  extensive  amount  of  work,  somehow  related  to  these 
concepts,  is  also  available  for  random  disorder  in  photonic  crystals  and  Anderson  localization  in 
random  materials  [1 1]-[16].  Here  we  concentrate  on  the  problem  of  periodic  arrays  of  plasmonic 
nanoparticles  in  order  to  derive  a  quantitative  expression  that  may  model  this  small  disorder  in 
closed  form,  with  a  first-order  perturbation  approach.  Different  from  the  other  approaches 
presented  in  the  literature  relative  to  this  problem,  our  analytical  results  allow  us  to  theoretically 
quantify  and  model  the  small  disorder  in  order  to  provide  simple  physical  analogies  and  insights 
that  could  be  later  applied  to  metamaterial  technology.  A  preliminary  portion  of  these  results  has 
been  presented  in  a  recent  conference  talk  [7]  and  these  results  have  been  recently  published  in 


c.  Theoretical  analysis 

Consider  an  infinite  linear  array  of  identical  particles  each  characterized  by  an  electric 
polarizability  aee  and  periodically  displaced  along  the  z  axis  at  distance  d  from  each  other 

(center-to-center  distance).  As  shown  in  various  recent  papers  on  the  problem  [17]-[23],  the 
guided  eigenmodes  supported  by  this  periodic  arrangement  propagate  along  the  array  with  an 
elp:  dependence,  where  [3  satisfies  the  following  equations: 


L : 

T : 


oo  _  _  _  _ 

6Z [ N~'d ~3  cos ( N(3d ) em  ( 1  -iNd) 

N= 1 

co  _  _ _  _ 

- 3^ I  AT3 d ~3  cos ( N/3 d ) em  ( 1  -iNd 

N=\  L 


(1) 
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where  d  =  k0d  ,  /?  =  J3  /  k0,  aee  =  klaee  /  ( 6  m:{) ) ,  k0=  oj^£nJult  ,  and  s0 ,  ju0  are  the  background 

medium’s  permittivity  and  permeability,  respectively.  In  Eq.  (1)  the  nanoparticles  fonning  the 
array  are  modeled  as  polarizable  dipoles.  The  two  equations  refer  to  the  two  orthogonal 
polarizations  for  the  eigenmodes  of  this  system,  i.e.,  respectively,  longitudinal  (Fig.  la)  and 
transverse  (Fig.  lb)  with  respect  to  the  axis  of  propagation  [9]. 

Eq.  (1)  supports  the  possibility  of  no-damping  sub-diffractive  guided  propagation  in  the  limit  of 
lossless  plasmonic  particles,  since,  as  expected,  the  mathematical  problem  represented  by  (1) 
ensures  a  real-valued  solution  for  any  \  <  p  <  n  I  d  under  the  lossless  condition  Im^a"1]  =  -1 

[9].  Therefore,  similar  to  any  periodic  metamaterial,  this  geometry  in  the  ideally  periodic  limit 
also  ensures  that,  despite  the  scattering  losses  from  each  one  of  the  nanoparticles,  the  overall 
behavior  of  the  chain  may  have  zero  radiation  losses. 

The  previous  equations  may  fully  take  into  account  the  general  case  in  which  the  particles 
forming  the  array  are  lossy.  In  such  case,  complex  solutions  for  p  are  required,  with  a 
corresponding  decay  in  the  direction  of  propagation.  The  summations  in  (1)  are  not  strictly 
convergent  in  this  situation  [19],  and  an  analytical  continuation  employing  polylogarithm 
functions  is  required,  obtaining  the  closed-fonn  dispersion  equations  first  introduced  in  [9]: 


L:3d 


=  aj 


3^-3 


T :  — d 
2 


f3(P,d)-id  f2(p,d) 
-f3(p,d)-idf2(p,d)-d2p(p,d) 


=  aj 


(2) 


where  d^  =  LiN[e'^+^d  ^  +  LiN{e  j  and  LiN  (x)  =  £  N  ^  \lt ,  with 

Lil  (x)  =  -ln(l-x)  is  the  polylogarithm  function,  as  defined  in  [24]. 


Fet  us  suppose  now  that  each  of  the  particles  in  the  linear  array  is  not  necessarily  positioned  at 
the  designated  location  zN  =  N-d ,  but  slightly  shifted  as  zN  =  Nd  +  SN ,  where  the  set  of  real 

numbers  SN  <sc  d  represents  quantitatively  an  imperfection  or  disorder  in  their  position  (see  Fig. 

1).  Analogously,  intrinsic  fabrication  limitations  may  generate  uncontrollable  small  differences 
among  the  particles,  and  each  of  the  polarizabilities  involved  in  the  summations  (1)  may  be 
rewritten  as  aN  =  aee  (1  +  cpN ) ,  with  cpN  being  a  set  of  complex  numbers. 


Po 
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Figure  1  -  Geometry  of  the  propagation  along  1-D  linear  chains  of  particles  under  the  two  orthogonal  polarizations 
[(a)  longitudinal,  (b)  transverse],  considering  random  imperfections/disorders  in  the  positioning  of  the  particles. 


In  order  to  model  the  presence  of  disorder  in  this  problem,  we  may  assume  that  these 
imperfections  are  inherently  random  in  nature  and  provide  small  perturbations  from  the  ideal 
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design  values.  Under  these  assumptions,  it  may  be  reasonably  expected  that  such  quasi-periodic 
system  may  still  be  able  to  guide  a  mode  with  e,ft:  space  dependence,  where  the  disorder  is 
anticipated  to  produce  an  analogously  small  perturbation  in  the  guided  wave  number  fj  .  Here  we 
aim  to  quantify  this  variation,  and  relate  it  to  the  expected  degree  of  disorder. 

The  dispersion  equations  (1)  have  been  obtained  for  the  two  polarizations  in  the  ideally-periodic 
geometry  by  assuming  that  the  induced  dipole  moments  along  the  chain  are  in  the  form 
p  v  =  p0£y/M  </ ,  with  N  being  any  integer  (positive  or  negative)  and  p0  being  the  dipole  moment 
induced  on  the  particle  located  at  the  origin.  Each  p^  produces  an  electric  field  at  the  origin, 
with  amplitude  E  N  =  G  ( Nd )  ■  p  v  ,  where  G  (  Nd )  is  the  dipolar  dyadic  Green  function.  Imposing 
that  the  sum  of  the  E  v  would  give  rise  to  a  total  electric  field  at  the  origin  self-sustaining  a 
dipole  moment  p0  provides  the  dispersion  equations  (l)-(2),  which  may  be  solved  for  the  guided 
wave  number  /3  once  the  array  geometry  is  known.  The  presence  of  small  disorder  in  the 
position  of  the  particles  has  a  perturbative  effect  to  the  original  dispersion  equation,  which  may 
be  analyzed  under  a  first-order  approximation  as  the  perturbation  to  the  self-sustaining  field  at 
the  location  of  p0 .  In  this  sense,  the  small  disorder  affects  the  evaluation  of  the  dispersion 

equations  in  two  distinct  ways:  the  Green  function  is  modified  into  G  ( Nd  +  8 N ) ,  due  to  the 

variation  of  the  distance  of  each  of  the  particle  from  the  origin,  which  causes  an  estimated 
variation  in  the  induced  field  on  the  particle;  moreover,  the  value  of  each  of  the  dipole  moments 
in  the  chain  is  also  varied  as  p^  =  ^0e^Nd+s^  by  the  variation  in  the  particle  position.  The 

variation  in  the  shape  or  electromagnetic  properties  of  the  particles  may  also  be  embedded  in  the 
estimation  of  the  induced  dipole  moment  of  each  particle,  which  corresponds  quantitatively  to 

the  expression  p  v  =  (l  +  ^)p0e'^Jva+<^ .  Summing  all  these  effects,  each  tenn  in  the  summations 

in  (1),  therefore,  takes  the  form: 

L  :  ( \N\  d +  SN  )~3  \  _  / 1 jy|  d  -  iSN )  ( 1  +  <pN ) 

(3) 

T :  ( \N\ d  +  SN y3  eWNl+s"' V(l^+^  ( 1  - / (|iV|  J +^)-(|iV|^  +  ^)2  j(1  +  ^)  ’ 

where  N  goes  from  — oo  to  go  and  it  takes  into  account  of  the  contribution  from  each  of  the 
particles  in  the  chain.  In  this  modified  form,  the  sums  in  the  dispersion  relations  cannot  be 
evaluated  in  any  closed  form,  since  each  of  their  terms  depend  separately  on  random  occurrences 
of  SN  and  <pN . 


Rearranging  the  terms  in  the  summation  and  expanding  in  Taylor  series  for  SN ,  owing  to  the 
assumption  of  small  disorder,  we  may  formally  rewrite  the  dispersion  relations  as  in  Eq.  (2): 


L:3d 


-3 


=  Vl 


3^-3 


T :  — d 
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fi{P’d)-id  f2{p,d ) 
f3(fi,d)-id  f2(]3,d)-d2f,(]3,d ) 


=  a  ' 


(4) 


where  the  following  expressions  for  the  modified  polarizabilities  hold: 
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(5) 


aL  ~  aee  ^  CN(Pn  +  +  +  /n^N"’  ' 

T  N=-cc 

The  coefficients  cN,  dN,  fN  are  the  Taylor  coefficients  of  (3).  As  an  example,  in  the 
longitudinal  polarization  they  take  the  form: 

cy=3(|JV|rf)"3e,We«(l-i|JV|rf) 

dN=-l>(Ndy  eip[Nd]e^d]  3-i[p  +  3) Nd  ~(/3  +  l)(Nd)2  .  (6) 

fN=3(Nd )~5  dp[m)e[Nd]  1 2 - 6i [fi  +  2) Nd  - (fi  + 1) (/?  +  5) ( Nd  f  +  i (/? + 1)2  ( Nd )3 

Eq.  (4)-(5)  combine  the  effects  of  disorder  along  the  chain  in  a  single  compact  expression, 
accurate  enough  to  evaluate  the  perturbation  of  the  dispersion  relations  due  to  the  random 
disorder  in  the  chain  in  terms  of  a  change  in  the  “average  polarizability”  of  the  particles  forming 
the  array,  allowing  a  quantitative  evaluation  of  the  effects  of  the  disorder  on  the  guided  wave 
number.  In  other  words,  the  disorder  along  the  chain  is  “seen”  by  the  propagating  mode  as  a 
variation  in  the  polarizability  of  the  particles  in  an  equivalent  ideally-periodic  array. 

Due  to  the  assumption  of  small  disorder,  the  summation  in  (5)  is  expected  to  weakly  perturb  the 
solution  of  Eq.  (4),  and  this  is  why  the  approach  of  summing  the  fields  at  the  particle  location  is 
still  applicable.  In  particular,  in  the  limit  of  interest  in  which  the  particles  are  lossless  or  with 
low-loss  (allowing  low-damping  modal  propagation  along  the  chain)  it  is  straightforward  to 
show  that  the  perturbation  in  (5)  principally  affects  the  imaginary  part  of  Eq.  (4),  which  is 
identically  satisfied  in  the  limit  of  lossless  particles.  This  implies  that  a  reasonably  small  disorder 
weakly  affects  the  phase  velocity  along  the  chain,  but  it  may  affect  more  considerably  its 
damping  coefficient  (imaginary  part  of  fd  ). 

Assuming  that  the  random  quantities  SN  and  <pN  are  Gaussian  distributions  with  zero  expected 
value  and  with  variances  a2s  and  cr ,  respectively,  we  may  evaluate  this  expected  variation  on 

the  effective  polarizability  of  the  particles  as  “seen”  by  the  mode.  The  Gaussian  hypothesis  on 
the  distribution  of  SN  and  <pN  represents  a  realistic  assumption,  since  the  disorder  associated 
with  technological  limitations  and  imperfections  may  be  reasonably  assumed  to  be  a  stochastic 
process  in  which  each  occurrence  of  SN  and  cp,  is  fairly  independent  from  the  others.  Due  to  the 

zero  mean  value  of  dN  and  cpN  and  their  independence,  the  dominant  term  in  the  summation  (5) 

for  evaluating  the  expected  values  and  (aT1^  is  given  by  fNSp .  The  expected  values  may 

be  evaluated  in  closed  form  using  the  variance  theorem  for  the  sum  of  Gaussian  distributions  and 
the  properties  of  the  polylogarithm  function.  Their  expressions  in  closed  form  are  given  in  the 
general  case  as: 
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3  j 

2  d5 


l2fs(fi,d)-6id(]3  +  2)LiA  e{p+1)d  )  +  6/ J  (/?  -  2)  Li  A  ei[M)d )  + 

-  d2  ( p  + 1 )  ( p  +  5 )  Z/3  ( el{P+')d  )  -  d 2  ( /?  - 1 )  ( /?  -  5 )  Z/3  (  e'i(p~l)d )  + 

+id 3  (/?  +  l)2  Li2  (eiP+l)d  )  +  /d"3  (/?  -  ^  Lii  [ei{P~')d ) 

12 f5  (jj,d)-6id(/3  +  2) Li4  iAp+x)d U  6 id  (p- 2)Li4  )  + 

-d2(p2+6]3  +  7)1/3  (el(p+')d )  - d2  ( ]32  - 6/3  +  l)L/  (e<M)d )  + 
+id 3  (A  + 1) []}  +  3 )z/2  (e'(/i+lH )  +  id 5  ( p  - 1) ( ^  -3 )Li2  (e'i{M)d )  + 
+d4(P  + 1)'  Lix  (g  W )  +  j4  ( A  - 1)'  Lix  [e^d ) 


(7) 


where  we  have  applied  the  properties  of  the  polylogarithms  to  evaluate  the  summation  (5)  in 
closed  form.  Even  if  the  fonn  of  these  expressions  is  quite  complex,  these  solutions  provide 
some  interesting  insights  into  the  effect  of  disorder  in  these  plasmonic  chain  waveguides. 

It  is  seen  in  (7)  that  the  main  perturbation  on  the  effective  polarizability  is  given  by  the  variance 
of  the  disorder  in  the  position  of  the  particles,  as  shown  by  (7).  Random  (and  independent) 
variations  in  the  geometry  of  the  particles  appear  to  compensate  each  other  in  the  first-order 
approximation  along  the  chain,  not  significantly  affecting  the  expected  value  (a  1 ) ,  as  it  was 

also  verified  numerically  in  [4]  for  3D  metamaterials  in  a  different  geometry.  The  variance  of  the 
position  disorder,  however,  which  is  a  measure  of  the  degree  of  disorder  in  the  system,  may 
significantly  perturb  the  effective  “average”  polarizability  of  the  particles  as  seen  by  the  guided 
mode,  since  the  scattering  losses  from  each  of  the  particles  cannot  be  completely  “canceled”  by 
the  other  particles  forming  an  ideally  periodic  array,  producing  residual  scattering  losses.  For  this 
reason,  in  the  following  sections  we  concentrate  on  the  analysis  of  random  disorder  in  the 
position  of  the  nanoparticles,  assuming  for  simplicity  that  they  have  the  same  size. 

Due  to  the  nature  of  Eq.  (4),  when  low-loss  particles  and  low-damping  guided  modes  are 
considered,  the  main  effect  in  the  perturbation  of  (7)  is  in  fact  expected  to  be  in  the  imaginary 
part  of  (7),  which  is  a  measure  of  the  “losses”  experienced  by  the  mode.  Similarly  to  what  was 
done  in  [9],  the  properties  of  the  polylogarithm  functions  allow  a  direct  evaluation  of  this 
imaginary  part.  Using  the  properties: 
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the  following  exact  relations  are  obtained  for  the  imaginary  parts  of  the  expected  values  of 
effective  polarizabilities  in  the  two  polarizations: 
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(9) 


under  the  assumption  that  Im[/?]  is  negligible  (low-damping  guided  modes).  This  result  is 

formally  consistent  with  the  first-order  approximation  for  the  effect  of  disorder  in  photonic 
crystal  waveguides  derived  in  [14]. 

We  reiterate  that  Eq.  (9)  applies  only  in  the  limit  of  small  disorder,  for  which  these  results 
represent  small  perturbations  of  the  solution  in  the  ideally  periodic  scenario.  In  this  limit,  this 
result  is  of  great  interest  for  its  inherent  simplicity.  In  this  context,  it  is  noticed  that  this 
perturbative  method  is  consistent  with  iterative  techniques  that  take  into  account  small  disorders 
in  periodic  arrays,  for  which  a  recent  application  to  periodic  metamaterials  has  been  presented  in 
[6].  Although  our  method  would  numerically  yield  consistent  results  for  any  type  of  random 
disorder  (with  zero  mean  value),  the  advantage  of  our  technique  when  applied  to  Gaussian 
disorder  is  that  it  provides  a  closed-fonn  solution  for  the  level  of  additional  losses  caused  by  the 
disorder,  directly  related  to  its  variance,  as  shown  in  Eq.  (9).  This  provides  relevant  physical 
insights  into  the  effect  of  small  disorder  in  periodic  structures,  as  we  further  discuss  in  the 
following. 

Without  resorting  to  any  further  approximation,  and  considering  the  fully  dynamic  interaction  in 
the  infinite  disordered  array  of  particles,  Eq.  (9)  indeed  states  that,  independent  of  the  average 
distance  between  the  particles,  the  disorder  in  the  positions  of  particles  affects  the  modal 
propagation  by  adding  an  effectively  additional  loss,  due  to  the  incoherent  scattering  radiation 
from  the  particles,  proportional  to  the  variance  of  the  disorder.  Similarly  to  the  closed-fonn 
solution  obtained  in  the  ideally-periodic  chain  for  the  imaginary  part  of  the  dispersion  relations 
(4),  which  ensured  power  conservation  along  the  chain  [9],  and  for  which  j  =  -1 ,  i.e.,  the 

particles  are  lossless  and  no-damping  propagation  is  obtained,  here  we  get  another  equally 
simple  generalized  expression  that  takes  into  account  the  presence  of  disorder  along  the  chain 
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and  adds  another  fonn  of  radiation  damping,  quantitatively  related  to  the  degree  of  disorder  in 
the  chain. 

We  note  that  the  sign  of  the  additional  tenns  in  (9)  is  strictly  negative,  ensuring  that  the  disorder 
increases  the  damping  factor  of  the  modes  whatever  is  the  geometry  of  the  chain  and  the  nature 
of  the  disorder.  For  passive  particles,  in  fact,  Im^a”1  J  =  - 1  -  ah]s ,  with  a~'ss  >0  [9].  Moreover, 

Eq.  (9)  is  consistent  with  the  physical  expectation  that  an  increase  in  the  variance  of  the  disorder 
would  correspond  to  an  increase  in  the  damping  factor  of  the  guided  modes,  and  that  a  tighter 
(and  slower)  propagation  along  the  chain  (increase  in  j3  )  makes  the  propagation  more  sensitive 
to  such  imperfections.  Another  interesting  feature  of  this  result  is  associated  with  the  fact  that  the 
transverse-polarization  configuration,  which  ensures  backward  propagation,  is  slightly  more 
sensitive  to  the  disorder  than  the  longitudinal-polarization  propagation.  This  is  consistent  with 
our  previous  findings  [9]  that  showed  how  this  polarization  has  a  slightly  lower  bandwidth  and 
more  sensitivity  to  losses  than  the  forward-wave  longitudinal  propagation. 

d.  Numerical  examples  and  validation 

Figure  2,  as  a  first  example,  shows  the  variation  of  Im[/?J  with  the  standard  deviation  <7S  of  the 

position  disorder  for  a  chain  of  silver  particles,  calculated  using  (9)  and  the  results  in  [9]. 
Longitudinal  polarization,  which  is  the  most  appealing  for  guiding  purposes,  is  considered  here 
and  in  the  following.  Moreover,  we  concentrate  on  the  positional  disorder,  which  has  been 
proven  in  the  previous  section  to  be  the  most  significant  mechanism  of  radiation  loss.  The  chain 
is  formed  by  silver  spherical  nanoparticles  with  average  radius  a  =  10  nm  and  center-to-center 
average  distance  d  =  22  nm .  We  have  evaluated  this  chain  geometry  considering  experimental 
values  of  the  silver  permittivity  available  in  the  literature  [25],  and  considering  the  inherent 
absorption  and  frequency  dispersion  (black  solid  line)  of  silver  at  optical  frequencies.  We  have 
also  added  the  red  dashed  curve,  which  neglects  the  absorption  in  the  particles,  in  order  to  isolate 
the  effect  of  disorder  on  the  damping  of  the  mode.  The  curves  in  this  case  have  been  evaluated  at 
the  free-space  wavelength  A0  =380 nm ,  which,  following  the  analysis  in  [9],  ensures  the  best 

guiding  properties  of  the  mode  as  the  highest  ratio  between  real  and  imaginary  parts  of  (3 .  At 
this  wavelength  the  pennittivity  of  silver  is  sAg  =  -3.8  +  /0.19  [25], 

As  predicted  by  the  previous  analysis,  increasing  the  disorder  along  the  chain  induces  a 
perturbation  of  ft ,  mainly  reflected  in  its  imaginary  part.  Small  deviations  from  the  ideal 
position  of  the  particles  do  not  considerably  perturb  the  guiding  properties  of  the  chain,  and  the 
absorption  is  usually  dominated  by  the  inherent  ohmic  absorption  of  metallic  particles.  Still,  the 
sensitivity  to  disorder  is  well  modeled  by  the  theoretical  considerations  in  this  geometry. 

Figure  3  shows  a  similar  case,  but  for  larger  silver  particles  and  correspondingly  larger 
separation  between  the  centers  of  the  particles,  highlighting  the  difference  in  the  effect  of 
disorder  with  a  change  in  the  geometry  of  the  chain.  In  this  case  a  =15 nm,  d=33nm  and 

/t0  =  390  nm  ,  for  which  eAg  =  -3 .9  +  i  0. 1 9 . 

It  is  evident  how  in  this  second  scenario  the  effect  of  disorder  is  somehow  lower.  This  is  due  to 
the  fact  that  in  this  second  scenario  the  particles  are  larger,  and  therefore  less  affected  by  the 
inherent  Ohmic  losses  and  effect  of  disorder,  since  the  electric  field  is  less  concentrated  in  their 
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volume.  Despite  the  mode  being  more  confined  around  the  chain  (higher  Re[/?] )  this  second 

example  ensures  more  robustness  to  Ohmic  losses  and  to  disorder,  consistent  with  its  inherent 
robustness  reported  in  [9]. 


Figure  2  -  Variation  of  imaginary  (a)  and  real  (b)  parts  of  (3  as  a  function  of  the  standard  deviation  of  disorder 
along  a  chain  of  silver  nanoparticles,  considering  (solid)  and  neglecting  (dashed)  Ohmic  absorption.  The  radius  of 
the  particles  is  a  =  10 nm  and  average  distance  between  any  two  neighboring  particles  (center  to  center)  is 
d  =  22  nm  at  the  free-space  wavelength  A0  =  380  nm  .  Here  and  in  the  following  the  longitudinal  polarization  is 

considered. 

Figure  4,  as  a  confirmation  of  the  validity  of  the  previous  theoretical  results,  reports  the 
numerical  simulations  of  a  linear  chain  of  particles  with  some  disorder  in  which  the  location  and 
the  size  of  each  particle  has  been  perturbed  from  the  original  periodic  geometry  of  Fig.  2.  In 
particular,  we  have  assumed  a  Gaussian  random  noise  in  the  position  of  the  particles  with 
standard  deviation  crs  =0.5  nm  and  a  deviation  from  the  ideal  radius  of  each  particle  a  =  10  nm 

equal  to  cra  =0.5  nm .  The  figure  reports  the  calculated  hn[/?]  for  the  ideal  chain  without 

disorder  and  for  the  chain  with  added  disorder,  both  evaluated  using  Eq.  (4),  corresponding  to  the 
black  solid  line  and  the  red  dashed  line,  respectively.  Moreover,  we  have  reported  the 
calculations  for  a  random  chain,  obtained  by  applying  (1),  (3)  and  truncating  the  summation  to 
the  first  20  elements  (that  is  evaluating  the  eigen-wave  number  considering  the  contributions 
from  the  nearest  40  nanoparticles  around  the  origin).  As  we  already  mentioned,  the  summation  in 
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(1)  is  not  strictly  convergent  when  complex  propagation  constants  are  considered,  due  to  the  fact 
that  the  mode  amplitude  grows  exponentially  in  one  direction  of  the  chain.  This  is  consistent 
with  any  complex  eigenmode  solution  (leaky-waves  or  lossy  geometries),  and  this  problem  is 
usually  solved  with  proper  analytical  continuation  of  the  series  in  the  complex  domain,  as  in  Eq. 

(4). 


b)  [  nm  I 

Figure  3  -  Similar  to  Fig.  2,  but  for  a  =  15  nm  and  average  distance  d  =  33  nm  at  the  free-space  wavelength 

A0  =  390 nm  . 

Since  now  the  summation  is  constituted  by  random  variables,  and  it  cannot  be  solved  in  closed- 
form,  we  have  considered  its  solution  obtained  after  proper  truncation,  which  is  a  good 
approximation  to  the  exact  analytical  continuation  derived  in  (2).  In  this  scenario,  where  each 
particle  has  a  random  deviation  in  position  and  size,  the  closed-form  analytical  solution  cannot 
be  derived,  but  the  summation  truncated  to  the  first  20  elements  provide  a  good  agreement  with 
our  theoretical  results  that  embed  the  effect  of  disorder  in  the  closed-form  expression  (4)-(5),  as 
evident  from  the  figure.  The  residual  oscillation  with  frequency,  evident  in  Fig.  4,  is  related  to 
the  truncation  effect,  as  shown  by  evaluating  Eq.  (l)-(3)  in  the  case  of  no  disorder  and  of 
increased  losses  using  Eq.  (4)  (black  and  red  dotted  lines,  respectively).  It  is  evident  that  the  line 
corresponding  to  the  random  chain  (dotted  blue  line)  agrees  very  well  with  the  truncated 
dispersion  relation  obtained  using  the  previous  theoretical  results.  We  have  verified  similar 
agreement  when  we  change  the  number  of  elements  considered  in  the  truncated  series. 

Figure  5  shows  the  full-wave  numerical  simulations,  performed  with  commercial 
electromagnetic  software  [26],  simulating  the  chain  geometry  of  Fig.  2  with  (Fig.  5a)  and 
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without  (Fig.  5b)  considering  the  disorder  in  the  position  of  the  particles,  in  the  amount  of 
<7g  =  2  nm .  These  simulations  refer  to  the  wavelength  A0  =  400  nm ,  for  which  the  silver 

permittivity  is  sA  -  -4.56  +  i 0.22  [25].  In  particular,  the  figure  reports  a  snapshot  of  the 

instantaneous  magnetic  field  distribution  orthogonal  to  the  plane  of  the  chain.  It  is  evident  from 
the  figure  how  the  guided  mode  along  the  chain,  although  supported  in  both  cases,  decays  faster 
due  to  the  presence  of  disorder,  consistent  with  the  results  of  the  previous  theory  and  with  the 
presence  of  an  effective  additional  loss  factor  in  the  particles  forming  the  chains.  In  the  figure  we 
have  also  reported  the  plot  of  the  longitudinal  electric  field  amplitude  sampled  on  a  line  parallel 
to  the  chain  axis  at  a  distance  of  80  nm  .  Despite  the  sharp  variations  due  to  the  granularity  of  the 
spheres,  a  clear  exponential  decay  is  visible  in  both  cases  in  the  figure,  which  is  highlighted  by 
the  dashed  lines  in  Fig.  5.  We  have  extracted  the  effective  Im[/?J  from  these  simulations,  which 
is  given  by: 


no  disorder:  Im[/?J  =  0.22 
disorder:  Im[/?]  =  0.35 


(10) 


respectively  corresponding  to  the  black  (darker)  and  red  (lighter)  lines.  As  expected,  these 
extracted  values  of  Im[/?]  are  consistent  with  the  values  predicted  by  Eq.  (9),  despite  the 
simplified  dipolar  approximation  of  our  analytical  model. 


Figure  4  -  Dispersion  of  ImT/?]  for  a  guided  mode  along  the  chain  of  Fig.  2  in  the  three  cases  of:  ideally  periodic 

chain  (black  solid  line),  disordered  chain  using  this  theory  (Eq.  (9),  red  dashed)  and  disordered  chain  using  a  random 
array  and  Eq.  (l)-(3)  truncated  to  Nmax  =  20  (blue  dotted).  The  disorder  in  the  position  of  the  particles  has  a 
standard  deviation  crs  =  0.5  nm  ,  whereas  the  disorder  in  the  size  of  the  particles  is  embedded  in  the  standard 
deviation  for  their  radius  (7a  =  0.5  nm  .  For  comparison,  here  we  have  also  reported  the  solution  of  Eqs.  (l)-(3) 
truncated  to  Nmax  =  20  for  an  ideally  periodic  chain  (black  dotted)  and  for  the  disordered  chain  using  Eq.  (9)  (red 

dotted). 
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b) 


Figure  5  -  Magnetic  field  distribution  orthogonal  to  the  plane  of  the  chain  for  the  geometry  of  Fig.  2  in  the  case  of: 

(a)  ideally  (periodically)  ordered  particles  and  (b)  with  random  disorder  in  the  position  of  the  particles  with 
<je  =  2  nm  .  The  plot  reports  the  electric  field  distribution  sampled  on  a  line  parallel  to  the  chain  axis,  at  a  distance  of 
80 nm  (solid  lines)  and  the  corresponding  exponential  decays,  consistent  with  our  theoretical  model. 

e.  Conclusions 

The  results  in  this  manuscript  confirm  that  the  disorder  in  plasmonic  waveguides  and 
metamaterials  may  give  rise  to  additional  scattering  losses  that  may  be  quantified  as  a  function  of 
the  degree  of  disorder  present  in  the  periodic  lattice.  The  geometry  analyzed  here,  a  linear  array 
of  plasmonic  nanoparticles,  has  been  selected  for  two  reasons:  the  possibility  of  a  closed-form 
analytical  solution,  generalizing  the  results  of  [9],  and  its  inherent  property  that  this  geometry 
shares  with  the  more  general  class  of  periodic  metamaterials,  i.e.,  supporting  low-loss 
propagation  by  the  means  of  cancellation  of  the  radiation  losses  that  an  ideally  periodic  lattice 
ensures.  It  has  been  proven  in  various  experimental  setups  that  sub-diffraction  optical 
waveguides  like  the  one  considered  here,  indeed  suffer  from  the  presence  of  losses,  and  signals 
cannot  travel  over  few  wavelengths  in  realistic  scenarios  with  subwavelength  cross-sectional 
size.  The  results  presented  here  allow  a  proper  quantification  of  the  effects  of  small  random 
disorder  on  the  expected  propagation  distance  and  on  the  further  increase  in  damping  due  to 
radiation  losses.  We  believe  this  may  be  useful  for  the  design  and  fabrication  of  these 
waveguides  and  of  optical  metamaterials  based  on  similar  concepts. 

With  similar  arguments,  these  results  may  be  extended  to  the  3D  scenario  of  closely  packed 
nanoparticles  forming  backward-wave  metamaterials  [10]  and  to  the  more  general  class  of 
periodic  metamaterials  and  photonic  crystals.  We  may  quantify  the  amount  of  expected 
scattering  losses  in  a  realistic  plasmonic  array  or  metamaterial  affected  by  disorder  associated 
with  technological  limitations  and  imperfections,  justifying  how  the  measured  losses  in  such 
setups  are  often  larger  than  those  expected  from  purely  theoretical  predictions  obtained  under  the 
assumptions  of  ideally-periodic  configurations.  It  is  interesting  to  underline  that  our  theoretical 
results  show  quantitatively  how  to  a  first-order  approximation  small  variations  in  the  shape  or 
electromagnetic  properties  of  the  inclusions  are  less  important  than  an  analogous  disorder  in  their 
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position,  for  what  concerns  the  overall  electromagnetic  behavior  of  the  metamaterial  setup  and  in 
particular  its  damping  and  absorption  properties.  To  conclude,  these  results  represent  an 
important  step  in  quantifying  the  radiation  and  scattering  losses  induced  by  disorder  in  plasmonic 
waveguides  and  periodic  metamaterials. 


4.  Waveguiding  Properties  of  Plasmonic  Nanoparticles  and  Voids 

a.  Summary 

In  this  section,  we  show  our  efforts  in  comparing  different  nanoscale  waveguiding  geometries 
involving  plasmonic  materials  for  sub-diffractive  propagation  at  optical  frequencies.  Deriving 
closed-form  analytical  fonnulas  to  analyze  the  different  structures,  we  show  how  the  presence  of 
a  plasmonic  background  may  produce  robust,  highly  confined  guided  wave  propagation  as 
compared  with  the  dual  setups  involving  plasmonic  particles  in  a  transparent  background. 
Advantages  and  disadvantages  of  different  scenarios  for  realizing  right-handed  and  left-handed 
propagation  in  one-dimensional  (ID)  and  two-dimensional  (2D)  waveguides  are  highlighted  and 
discussed. 


b.  Introduction 

Plasmonic  phenomena  and  their  related  applications  have  witnessed  fantastic  development  in 
recent  years.  In  particular,  at  optical  frequencies  where  plasmonic  materials  are  readily  available 
in  nature  [27]-[28],  various  potential  applications  of  plasmonic  concepts  have  been  recently 
proposed,  including  transport  of  optical  energy  with  tight  lateral  confinement  [29]-[56],  left- 
handed  wave  propagation  [5 7] -[60],  and  slow-light  devices  [61]-[63],  just  to  name  a  few.  As  one 
such  application,  inspired  by  our  idea  of  bringing  the  notion  of  circuit  elements  into  optical 
frequencies  [64],  we  have  proposed  and  analyzed  different  geometries  for  realizing  optical 
nanotransmission-lines  in  the  form  of  chains  of  plasmonic  nanoparticles  [39]  and  nanorods  [41] 
(ID  propagation),  planar  plasmonic  nanolayers  [56]  (2D),  and  plasmonic  3D  nanoarrays  of 
particles  to  form  optical  nanomaterials  [40].  Bringing  the  notion  of  transmission-line  from  the 
radio  frequencies  (RF)  and  microwaves  into  optical  frequencies  may  inspire  the  design  of  guided 
modes  more  tolerant  to  Ohmic  absorption,  less  radiation  losses  and  larger  bandwidth  of 
operation,  all  inherent  properties  of  RF  transmission  lines.  The  different  waveguiding  geometries 
that  we  have  recently  proposed,  all  based  on  the  peculiar  properties  of  plasmonic  materials  at  IR 
and  optical  frequencies,  indeed  allow,  under  proper  conditions,  reasonably  long  propagation 
distances  with  highly  sub-diffractive  lateral  confinement  of  the  guided  beams  and  relatively 
larger  bandwidth  of  operation. 

Even  if  ideally  the  plasmonic  nanowaveguides  may  be  designed  to  totally  suppress  unwanted 
radiation  [3 9] -[40],  the  presence  of  small  disorder,  Ohmic  absorption  and  other  unwanted 
imperfections  may  generate  radiation  losses  in  the  transparent  background,  which  affect  the 
propagation  distance  of  the  guided  beams  [65].  For  this  reason,  we  have  recently  been  interested 
in  analyzing  in  detail  the  dual  of  these  setups,  i.e.,  analogous  geometries  employing  dielectric 
waveguides  surrounded  by  a  plasmonic  background  material,  in  the  form  of  arrays  of  voids  or 
dielectrics,  cylindrical  voids,  or  dielectric  nanolayers  embedded  in  noble  metals  or  other 
plasmonic  materials  as  backgrounds,  consistent  with  some  recent  analogous  proposals  [42]-[54]. 
Due  to  the  opaqueness  of  such  materials  at  optical  frequencies,  wave  propagation  in  the 
background  is  forbidden  and  the  corresponding  radiation  losses  of  the  guided  modes  may  be 
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reduced.  The  confinement  of  the  guided  beams  inside  the  dielectrics  or  voids,  which  usually  have 
lower  loss,  may  also  ensure  a  lower  sensitivity  to  Ohmic  absorption,  which  is  another  possible 
advantage  of  these  dual  configurations.  This  is  of  course  valid  as  long  as  comparable  field 
confinement  is  achieved,  since  large  field  concentration  would  always  lead  to  larger  sensitivity  to 
absorption  loss  in  both  setups. 

In  general,  a  common  trend  in  all  the  waveguiding  geometries  involving  plasmonic  effects  may 
be  emphasized:  a  trade-off  between  propagation  length  of  the  guided  waves  and  the  limitation  on 
its  lateral  confinement  and  bandwidth.  In  this  sense,  a  figure  of  merit,  defined  as  the  ratio  of 
propagation  distance  versus  beam  cross  section  may  be  defined  in  order  to  sort  out  the  most 
promising  waveguiding  geometries  involving  plasmonic  materials  at  the  frequency  of  interest. 
Moreover,  the  fact  that  the  power  inside  plasmonic  materials  generally  flows  in  a  direction 
opposite  to  the  modal  phase  velocity  [56]  allows  one  tailoring  these  waveguides  to  exhibit  either 
forward-  or  backward- wave  propagation  (i.e.,  the  group  velocity  being  parallel  or  anti-parallel 
with  the  phase  velocity  of  the  mode),  depending  on  whether  the  power  is  more  concentrated  in 
the  plasmonic  or  in  the  dielectric  region.  Efficient  and  optimized  backward-wave  waveguides 
may  play  an  important  role  in  the  framework  of  left-handed  material  research  and  sub¬ 
wavelength  imaging,  as  discussed  in  Ref.  [56]. 

In  this  section,  after  discussing  the  analytical  properties  of  the  guided  modes  in  these  different 
types  of  waveguides,  we  compare  and  contrast  their  properties,  exploring  which  ones  of  the 
different  proposed  waveguiding  geometries  and  solutions  may  be  more  robust  to  material 
absorption  and  radiation  losses,  and  which  ones  exhibit  a  “better”  bandwidth  of  operation. 


c.  Plasmonic  and  dielectric  cavities 


It  is  well  known  that  a  subwavelength  plasmonic  nanoparticle  of  radius  a  and  pennittivity  s 
may  support  a  “quasi-static”  resonance  when  Re[<?]<0  [28].  This  phenomenon  may  be 

interpreted  in  terms  of  the  nanocircuit  concepts  [64]  as  the  intrinsic  resonance  between  the 
nanoinductor  elements  represented  by  a  plasmonic  nanoparticle  and  the  nanocapacitor 
constituted  by  its  fringing  fields  in  a  background  material  with  pennittivity  sb  >  0 .  The  specific 

values  of  the  two  nanocircuit  elements,  as  described  in  Ref.  [64],  depend  on  the  material 
parameters  and  the  geometry  of  the  nanoparticle,  which  determines  the  combined  resonance 
condition.  It  is  well  known,  for  instance,  that  the  resonance  of  a  plasmonic  nanosphere  is 
obtained  when: 


=  0, 


(ID 


j«  (ka)  y„  {Ka) 

\ka  jn  (ka)]  Is  [kba  yn  (kbaj\  / , 

where  jn  and  yn  are  spherical  Bessel  functions,  depending  on  the  resonant  order  n  ,  k  -  co^£ju0 
,  kb  -  (with  the  choice  for  the  branch  root  of  kh  satisfying  the  radiation  condition  at 

infinity)  and  ju0  is  the  free-space  permeability  assumed  to  be  the  same  in  all  the  (non-magnetic) 
materials  in  our  study.  In  the  small-radii  approximation,  Eq.  (11)  is  satisfied  near  the  frequency 

Yl  +  1 

for  which  Re[^]  - - sh .  Different  resonant  conditions  may  be  achieved  for  different  shapes 
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of  the  nanoparticle  or  with  a  multi-layered  geometry  [67].  For  the  dominant  resonant  mode,  with 
n  =  1 ,  the  corresponding  resonance  Q-factor  may  be  calculated  as: 


0 


3 

2(V0'3+3^ 

Jo 


(12) 


where  we  have  assumed  a  Drude-model  dispersion  for  the  nanosphere  permittivity  as 


£  =  £n 


1- 


3/02 


f(f+ir ) 


with  f0  being  the  resonance  frequency.  (Here,  we  are  assuming  that  the 


plasma  radian  frequency  for  the  spherical  particle  is  co  -  2x^3 fa .)  The  denominator  in  Eq. 


(12),  proportional  to  the  total  losses  in  the  system,  is  evidently  given  by  the  sum  of  two  terms: 
the  radiation  losses,  inversely  proportional  to  the  nanosphere  volume,  consistent  with  Ref.  [66], 
and  the  ohmic  absorption  in  the  material,  proportional  to  the  damping  frequency  y  of  the 
material.  Even  in  an  ideally  lossless  particle,  the  Q-factor  in  Eq.  (12)  would  be  indeed  limited  by 
the  intrinsic  radiation  losses  of  the  particle.  For  higher-order  resonances  the  tenn  in  the 

denominator  of  (12)  becomes  ( k0a )  2"  ' ,  implying  smaller  radiation  losses  and  higher  Q,  as 


expected. 


If  we  now  reverse  the  role  of  inductors  and  capacitors,  i.e.,  we  analyze  the  dual  geometry 
composed  of  a  dielectric  nanosphere  with  pennittivity  e  >  0  embedded  in  an  s-negative  (ENG) 
background  with  permittivity  sb<  0 ,  the  resonance  is  described  by  an  analogous  equation: 


h  (ka) 

(M 

\kcijn{ka)\  Is 

kbah^(kba) 

!eb 

(13) 


where  is  the  spherical  Hankel  function  of  first  kind,  which  is  purely  imaginary  for  imaginary 
kb  (associated  with  the  fact  that  the  ENG  background  does  not  support  wave  propagation),  and 

the  field  distribution  is  mainly  concentrated  in  the  dielectric  nanosphere  and  at  the  plasmonic 
interface,  decaying  much  more  rapidly  in  the  outside  region.  In  this  case,  the  expression  for  the 
Q-factor  simplifies  to: 

0  =  — ,  (14) 

r 


f 


where  we  have  assumed  sb=s 


1- 


2n  +  1 
n  + 1 


A 


f2 

J  o 


,  which  once  again  supports  a  resonance  at 


V  J 

frequency  f0 .  Eq.  (14)  is  valid  for  any  resonant  order  n  . 


It  is  interesting  to  see  that  such  a  cavity  may  support  an  extremely  high-  Q  resonance,  limited 
only  by  the  inherent  losses  of  the  background  material,  since  no  radiation  is  involved  in  the 
unbounded  plasmonic  background.  Interestingly,  the  cavity  size  is  not  limited  by  diffraction,  i.e., 
in  principle  its  dimension  may  be  made  arbitrarily  small,  since  the  field  distribution  would  be 
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adjusted  in  order  to  remain  concentrated  around  the  plasmonic  interface,  with  the  Q  as  derived 
in  Eq.  (4),  independent  of  the  cavity  volume  or  of  the  resonant  order. 

Figure  6  reports  the  variation  of  the  electric  field  distribution  induced  on  the  surface  of  a  glass 
spherical  cavity  embedded  in  the  silver  background,  with  radius  a  =  5  nm ,  for  a  sample 
excitation.  It  is  interesting  to  note  how  the  different  resonant  peaks  1  have  similar  operating 
bandwidths,  since  they  are  characterized  by  a  similar  Q-factor,  as  predicted  by  Eq.  (14)  . 


Frequency  [THz  ] 

Figure  6  -  Electric  field  amplitude  on  the  surface  of  a  glass  cavity  with  radius  a  =  5  nm  embedded  in  a  silver 

background  (considering  realistic  losses  of  silver). 


d.  Chains  of  dielectric  voids  in  an  ENG  background 

If  the  0-D  problem  of  a  sub-wavelength  dielectric  cavity  in  an  ENG  background  presents  some 
interesting  peculiarities,  collections  of  such  small  cavities  may  give  rise  to  other,  even  more 
exciting  possibilities  [42],  [47],  [50].  For  instance,  consider  cascading  several  identical  dielectric 
voids,  of  any  arbitrary  shape  and  material  properties,  in  a  linear  array  format,  with  center-to- 
center  spacing  d  >2a  .  The  propagation  along  these  arrays  may  be  modeled  analogously  to  the 
previous  section  and  to  Ref.  [39],  where  we  have  solved  analytically  the  propagation  along  linear 
arrays  of  plasmonic  particles  in  free-space.  In  particular,  when  our  technique  in  Ref.  [39]  is 
applied  to  this  dual  scenario,  the  dispersion  equations  for  the  guided  modes  supported  by  such 
arrays  of  nanovoids  embedded  in  an  ENG  background  become,  in  the  limit  of  negligible  losses: 
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(15) 


1  As  seen  in  Fig.  1,  the  n  =  1  resonance  arises  around  7 55  THz  ,  the  other  resonances  arise  in  order  when  the 
frequency  is  lowered,  consistent  with  Eq.  (11). 

2  This  is  different  from  the  behavior  of  a  plasmonic  nanoparticle  embedded  in  a  dielectric  (e.g.,  air),  where,  due  to 
radiation  losses,  the  bandwidth  of  higher-order  resonances  is  usually  much  narrower  than  the  first-order  resonance. 
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Respectively,  for  the  longitudinal  and  transverse  polarization  (with  respect  to  the  array  axis  x , 
see  the  insets  in  Fig.  7a  and  Fig.  7b  respectively),  where  [3  is  the  longitudinal  wave  number  of 
the  guided  modes  ( e,ix  spatial  dependence)  and  aee  is  the  generic  electric  polarizability  of  each 

individual  nanovoid,  consistent  with  Ref.  [39].  Since  the  wave  propagation  is  forbidden  in  the 
background  material,  the  radiation  from  this  array  is  automatically  avoided  in  this  geometry,  and 
the  coupling  among  distant  elements  is  consistently  minimized.  In  this  lossless  limit,  Eq.  (15) 
indeed  ensures  quick  convergence  and  only  a  few  terms  should  be  considered  for  an  accurate 
prediction  of  the  guided  wave  number  f3  .  Introducing  losses,  however,  Eq.  (15)  is  not  properly 
converging  and  an  analytic  continuation  is  necessary  to  handle  the  problem  accurately.  In  this 
sense,  the  closed-form  expressions  may  be  derived,  analogous  to  those  derived  in  the  previous 
section  for  the  present  dual  geometry,  as: 
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with  fN(/3,d)  =  LiN{e<M^)  +  LiN{e-^),  LiN(z)  =  fj 0  =  0  /  (~ikb) ,  d  =-ikhd , 
-ik  3 

aee  = - —ocee.  LiN  is  the  polylogarithm  function,  which  may  be  analytically  continued  in  the 

6  neh 

complex  plane  following  Refs.  [39],  [68]. 

Owing  to  the  absence  of  radiation,  these  equations  are  much  simpler  to  analyze  than  their 
counterparts  derived  in  Ref.  [39]  for  the  transparent  background  scenario.  When  losses  are 
absent,  they  are  real-valued.  When  the  losses  are  considered  in  general  their  imaginary  part  is 
simply  associated  with  the  ohmic  absorption  in  the  background  (notice  that  due  to  the 
nonnalizations,  when  losses  are  considered  the  normalized  distance  d  becomes  a  complex 
quantity).  Due  to  the  properties  of  LiN  (z)  and  their  derivatives  (see  Ref.  [39]  for  details),  Eq. 

(16)  implies  that  in  the  longitudinal  polarization  the  chain  supports  only  “backward”  modes 
(oppositely-signed  phase  and  group  velocities),  whereas  the  modes  are  inherently  “forward”  for 
the  transverse  polarization.  This  is  opposite  to  the  dual  scenario  of  arrays  of  plasmonic  particles 
in  a  transparent  background,  consistent  with  the  fact  that  in  both  polarizations  we  have 
effectively  interchanged  the  roles  of  nanoinductors  and  nanocapacitors  (in  the  longitudinal 
polarization,  the  equivalent  circuit  model  for  this  geometry  consists  of  series  capacitors  and 
shunt  inductors,  like  in  a  left-handed  transmission  line[69],  whereas  for  transverse  modes  we 
have  series  inductors  and  shunt  capacitors,  as  in  a  regular  transmission  line  supporting  forward 
wave  propagation).  Eq.  (16)  detennines  the  conditions  on  the  polarizability  of  each  individual 
cavity  and  on  the  spacing  among  them,  so  that  such  arrays  may  support  propagation.  In 
particular,  the  range  of  guidance  for  the  two  polarizations  may  be  written  in  closed-form  by 
considering  the  properties  of  LiN  (z)  and  the  fact  that,  due  to  its  periodicity,  the  array  may 

support  propagation  in  the  range  0  <  (3  <  n  .  The  range  of  guided  wave  propagation  as  a  function 
of  aee  may  be  interestingly  written  in  closed  form  as: 
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valid  in  the  limit  of  zero  losses. 
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Figure  7  -  Regions  of  guidance  (between  the  black  and  red  lines)  for  arrays  of  nanovoids  in  an  unbounded  plasmonic 
background  for  the  two  polarizations,  as  a  function  of  the  normalized  distance  d  for  fixed  ratio  d  /  a  =  2.1. 

These  results  are  valid  for  any  shape  of  the  nanovoids,  since  they  are  written  in  terms  of  the 
generic  polarizability  aee.  For  a  homogenous  spherical  void  of  permittivity  s,  we  have 

a~xe  -  — ikba )  3  £  +  ^£b  and  this  allows  plotting  the  conditions  (17)  in  terms  of  the  ratio  sb!  e , 


as  reported  in  Figure  7  for  the  case  of  d  =  2.1a  .  The  solid  lines  in  the  figure  are  relative  to  the 
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current  scenario  (i.e.,  voids  embedded  in  the  ENG  background),  whereas  the  dashed  lines  refer  to 
the  dual  case  of  plasmonic  particles  in  a  transparent  background  (consistent  with  the  analogous 
plots  in  Ref.  [39]).  The  regions  of  guidance  are  delimited  by  the  loci  /?  =  0  and  /?  =  n  /  d 
(notice  that  in  this  geometry,  for  which  propagation  is  forbidden  in  the  background  material,  also 
fast  modes  with  /?  <  1  may  propagate  without  losses,  like  inside  a  regular  waveguide  at  RF, 
differently  from  the  case  of  plasmonic  particles  that  are  required  to  radiate  away  when  J3  <kh). 
In  the  plot  other  curves  for  specific  values  of  fd  within  the  guidance  region  are  also  shown. 
Since  £h<0  and  dsht dco>  0  for  passivity  requirements  [70],  both  plots  confirm  that 

longitudinal  modes  are  necessarily  backward  (Fig.  7a)  and  transverse  modes  are  forward  (Fig. 
7b),  opposite  to  the  dual  configuration  in  Ref.  [13]. 

Some  other  interesting  features  are  evident  in  these  plots:  in  principle  there  is  no  limit  on  the 
distance  between  the  nanovoid  cavities  to  support  propagation,  even  if,  for  large  enough 
distance,  propagation  is  allowed  only  very  near  the  resonant  frequency  of  the  individual  cavities 
(for  which  sb  =  —s  /  2  ).  This  may  substantially  enlarge  the  relative  region  of  propagation  when 

compared  with  the  dual  configuration  of  plasmonic  particles  (dashed  lines),  particularly  for  the 
transverse  polarization  (in  Ref.  [39]  it  was  proven  analytically  that  arrays  of  plasmonic  particles 
in  a  transparent  background  may  support  propagating  modes  with  no  radiation  losses  only  for 
d  <  n  ).  It  should  be  noticed,  however,  that  for  large  distances  the  supported  modes  are  very  fast 
( fd  —  0 )  and  this  propagation  region  may  become  sensitive  to  the  material  absorption  in  the 
background.  For  closer  and  closer  particles,  the  range  of  values  of  permittivities  of  dielectric 
voids  for  which  propagation  is  admissible  gets  relatively  larger,  consistent  with  the  tight 
coupling  among  the  individual  resonances  in  the  chain,  and  therefore  relatively  larger 
bandwidths  of  propagation  are  expected,  compared  to  the  dual  geometry,  due  to  absence  of 
radiation  loss  from  the  array.  In  the  limit  of  d  — »  0  ,  the  regions  of  guidance  coincide  for  the  two 
dual  geometries,  since  in  this  case  quasi-static  considerations  hold  [39]. 

Figure  8,  as  an  example,  shows  the  dispersion  of  the  guided  modes  for  a  chain  of  spherical  SiC 
voids  embedded  in  a  silver  background  with  a  =  5nm  and  d  =  10.5 nrn ,  considering  realistic 
values  of  pennittivity  (including  dispersion  and  losses)  for  both  materials  [27], [71].  In  particular, 
here  the  silver  is  modeled  using  a  Drude  model  with  plasma  frequency  /  =2. 175 Fife  and 

damping  frequency  y  =  43 5 Fife.  The  figure  reports  the  real  and  imaginary  parts  of  the 
nonnalized  guided  wave  number  Jd  /  k0,  with  k0  =  co^soJu0  .  Fig.  8a  reports  the  comparison 

between  the  dispersion  curves  in  the  ideally  lossless  scenario  (neglecting  loss  in  silver,  solid 
lines)  and  the  case  in  which  losses  are  considered  (dashed).  As  expected,  their  influence  is  most 
visible  near  the  edges  of  the  pass-band,  since  the  group  velocity  is  lower.  The  thin  dashed  line  in 
this  panel  represents  the  dispersion  of  the  curve  Re[/?]  =  n  Id  ,  which  coincides  with  the  stop- 

band  for  these  guided  modes,  as  again  expected.  Fig.  8b  reports  the  imaginary  part  of  the 
nonnalized  guided  wave  number  in  the  two  polarizations,  also  reporting  the  level  of  absorption 
in  dB  I  / 10,  where  A0  =  In  /  k0 .  It  is  clear  that  realistically  these  guided  beams  cannot  propagate 
over  multiple  wavelengths,  but  this  is  expected,  considering  the  level  of  confinement  represented 
by  the  large  value  of  Re[/?] . 
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Figure  8  -  Dispersion  of  the  real  (a)  and  imaginary  (b)  part  of  P  ,  normalized  to  k0 ,  for  a  linear  array  of  SiC  voids 
of  radius  a  =  5  nm  in  a  silver  background  with  d  =  10.5  nm  .  The  thin  dotted  line  corresponds  to  the  condition 
P  =  n  /  d  ,  which  is  the  upper  limit  for  modal  propagation,  whereas  the  dashed  lines  refer  to  the  case  in  which 

material  losses  are  neglected. 

It  is  evident  that  the  dispersion  of  p  versus  frequency  supports  the  backward  (forward) 
propagation  of  longitudinal  (transverse)  modes,  with  the  derivative  dp  /  dco  confirming  the 
previous  discussion.  The  bandwidth  for  longitudinal  modes  is  slightly  larger  than  for  the 
transverse  modes,  consistent  with  the  dual  scenario  [39],  implying  that  in  this  geometry 
backward-wave  propagation  is  actually  more  robust  than  forward-wave  propagation,  in  terms  of 
bandwidth  and  losses  (Fig.  8b).  The  possibility  of  relatively  low-loss  propagation  over  a 
relatively  wide  bandwidth  is  feasible,  consistent  with  the  transmission-line  analogy  at  radio- 
frequencies.  Another  advantage  of  this  dual  configuration,  with  respect  to  plasmonic  particles  in 
a  transparent  background,  is  the  suppression  of  the  spurious  transverse  mode  with  forward 
properties  [39]  that  limits  the  bandwidth  of  the  transverse  backward- wave  mode  in  the  dual 
geometry.  Here  its  presence  is  avoided  due  to  the  absence  of  propagation  in  the  background 
material.  In  the  figure  the  dashed  lines  refer  to  the  case  in  which  material  losses  are  neglected, 
exhibit  the  fact  that  their  presence  does  not  sensibly  affect  the  phase  velocity  of  guided  modes 
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(Re[/?] ,  Fig.  8a),  but  it  simply  affects  the  propagation  distance  (Im[/?] ,  Fig.  8b).  As  noticed  in 

panel  a),  and  consistent  with  the  dual  scenario  [39],  losses  are  more  significant  in  the  regions 
where  slow-wave  propagation  arises,  i.e.,  near  the  edges  of  the  propagation  band,  where  the 
group  velocity  deal dfd  -  0  . 


e.  Dielectric  nanorod  in  an  ENG  background 

In  the  limit  in  which  the  dielectric  nanovoids  analyzed  in  the  previous  section  are  closely  packed 
together,  they  may  constitute  a  homogeneous  ID  cylindrical  nanorod  of  radius  a  surrounded  by 
an  ENG  background,  as  investigated,  e.g.,  in  Ref.  [55].  The  properties  of  the  dual  configuration, 
i.e.,  a  plasmonic  nanorod  used  as  waveguide,  have  been  thoroughly  investigated  by  various 
groups  in  several  papers  (e.g.,  Refs.  [41], [72])  and  the  analytical  framework  is  the  same  when 
plasmonic  materials  form  the  background.  In  particular,  the  dispersion  equation  for  a  mode  with 
generic  azimuthal  order  n  supported  by  a  nanorod  with  pennittivity  s  and  radius  a  is  given  by: 
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Still  sub-diffraction  propagation  may  be  obtained  and,  not  surprisingly,  the  properties  of  the 
chain  configuration  analyzed  in  the  previous  section  do  appear  in  the  modes  guided  by  these 
cylindrical  nanovoids.  In  particular,  the  longitudinal  polarization  in  the  array  of  voids 
corresponds  to  the  azimuthally  symmetric  TM0  mode  (both  with  backward-wave  properties), 

whereas  the  transverse  polarization  in  the  chain  corresponds  to  the  forward  hybrid  quasi- EM, 

mode,  consistent  with  the  results  in  Ref.  [41].  The  important  difference  between  these  two 
geometries  arises  in  the  frequency  band  in  which  propagation  is  supported.  If  the  array  of 
nanovoids  supports  propagation  in  a  frequency  band  centered  around  the  resonance  of  each 
nanovoid  (for  the  homogeneous  sphere  case  when  eb  —  —s  12),  the  local  geometry  is  drastically 

changed  when  the  array  of  inclusions  merges  into  a  single  infinite  nanorod,  shifting  the 
propagation  band  around  the  plasmonic  frequency  of  the  nanorod,  which  happens  for  a 
homogeneous  cylinder  near  sb  =  —s  for  any  order  n  >  0 .  This  is  shown  in  Fig.  9,  where  we  have 

reported  the  propagation  properties,  in  both  polarizations,  of  a  cylindrical  nanorod  made  of  SiC 
in  a  silver  background,  with  a  =  5 nm  ,  consistent  with  the  example  in  Fig.  8.  The  plots  confirm 
the  correspondence  between  the  cylindrical  guided  modes  and  polarization  of  the  modes  in  the 
arrays  described  in  the  previous  section,  showing  how  the  cylindrical  geometry  may  provide 
another  possibility  for  obtaining  highly  confined  guided  wave  propagation  in  specific  frequency 
bands  of  interest  with  relatively  long  propagation  distances. 
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Figure  9  -  Dispersion  of  the  real  (a)  and  imaginary  (b)  part  of  P  ,  normalized  to  Re[-/k6]  calculated  at  the  central 
frequency,  for  a  cylindrical  nanorod  made  of  SiC  with  radius  a  =  5  nm  embedded  in  an  unbounded  silver 

background. 

We  compare  more  thoroughly  the  guidance  properties  of  the  ID  waveguides  analyzed  in  these 
two  subsections  in  the  following.  It  would  be  interesting  to  compare  more  thoroughly  the 
propagation  along  nanovoid  arrays  and  the  dual  geometry  of  plasmonic  nanoparticle  arrays,  as 
well  as  nanorods  and  cylindrical  nanovoids.  In  part,  we  also  address  this  comparison  in  the 
following  section.  However,  it  is  worth  noting  that  it  is  not  straightforward  to  establish  a  fair  and 
universal  criterion  to  compare  these  different  waveguide  geometries.  On  the  one  hand,  a 
thorough  comparison  would  require  keeping  the  same  shape  of  the  waveguide  and  using  the 
same  materials,  with  similar  level  of  losses.  However,  as  it  is  clear  after  comparing  the 
dispersion  bands  in  Fig.  7,  the  guidance  regions  in  these  two  complementary  scenarios  would  fall 
in  very  different  frequency  bands.  For  instance,  for  the  longitudinal  mode  of  the  arrays  of 
nanovoids  in  Fig.  8,  propagation  is  expected  in  the  green-blue  bands,  whereas  the  dual  setup  of 
silver  nanoparticles  in  a  SiC  background  would  support  propagation  in  the  red  spectrum.  This 
drastic  difference  is  of  course  reflected  in  substantially  different  properties  of  the  involved 
plasmonic  materials,  which  affect  the  overall  propagation  lengths.  Similar  considerations  are 
applicable  in  the  case  of  cylindrical  rods. 

f  2D  arrays  of  nanovoids 

The  extension  of  the  analytical  solution  presented  in  Ref.  [39]  and  above  to  the  2D  array  of 
nanoparticles  and  nanovoids  may  be  performed  following  the  analytical  technique  that  we  have 
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recently  applied  to  study  the  propagation  in  3D  arrays  of  plasmonic  particles  [40],  based  on  the 
acceleration  technique  presented  in  Ref.  [73].  In  this  2D  geometry,  Eq.  (16)  becomes: 
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where  J3x,j3y  are  the  components  of  the  propagation  wave  vector  in  the  x-y  array  plane  and 
dx ,  d  are  the  spacing  distances  (i.e.  the  periodicities  in  the  x  and  y  directions).  These  fonnulas 
are  equally  valid  both  for  2D  arrays  of  plasmonic  nanoparticles  and  of  nanovoids  in  a  plasmonic 
background,  once  again  in  terms  of  the  individual  polarizability  of  each  inclusion  a~l .  The 
convergence  of  the  summations  in  (19)  is  very  fast  (few  terms  are  needed  for  convergence)  and 
the  equations  are  also  valid  for  complex  values  of  /?  (due  to  ohmic  absorption  or  radiation  losses 
in  the  transparent  background  case). 

The  properties  of  these  chains  are  analogous  to  those  highlighted  in  Section  2b  for  the  ID  case, 
even  if  the  operational  bandwidth  and  robustness  to  losses  may  be  worsened  in  the  transverse 
polarization  when  the  parallel  chains  get  too  close  to  each  other,  consistent  with  our  findings  in 
3D  arrays  [40],  In  the  case  of  nanovoids  in  a  background  material,  in  this  2D  case  the  transverse 
propagation  also  implies  forward-wave  modes  and  longitudinal  propagation  corresponds  to 
backward  propagation.  Longitudinal  modes  are  generally  more  robust  to  losses  and  have  a 
slightly  larger  frequency  bandwidth,  and  in  both  scenarios  the  propagation  band  is  centered 
around  the  resonant  frequency  of  the  individual  inclusions,  similar  to  the  ID  arrays. 
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g.  Planar  nanolayers 


Similar  to  the  ID  geometry,  when  the  nanoparticles  or  nanovoids  are  closely  packed  together,  the 
2D  structure  resembles,  in  the  limit,  a  single  planar  nanolayer  with  thickness  2 a ,  corresponding 
to  the  geometry  thoroughly  analyzed  in  Ref.  [56].  The  corresponding  dispersion  relation  for  the 
TM  modes  may  be  split  into  even  and  odd  operation,  with  conditions: 
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As  we  have  shown  in  Ref.  [56],  this  geometry  also  fully  supports  the  nanocircuit  interpretation, 
acting  as  an  infinite  cascade  of  nanoinductors  and  nanocapacitors  arranged  in  series  or  in  parallel 
as  a  function  of  the  polarization  of  the  local  electric  field.  In  particular,  the  longitudinal 
polarization  in  the  2D  nanochain/nanovoids  problem  corresponds  to  the  odd  TM  modes  in  the 
planar  layer  geometry,  and  correspondingly  transverse  polarization  corresponds  to  the  even 
modes.  This  is  confirmed  by  the  backward  (forward)  nature  of  the  odd  (even)  modes  in  the 
metal-insulator-metal  geometry  [56],  Moreover,  in  this  case  the  propagation  band  is  generally 
located  around  the  resonance  of  the  layer  (when  £b  =-e),  even  if  usually  the  bandwidths  over 
which  propagation  takes  place  in  this  planar  layer  geometry  are  very  large.  In  the  following 
section,  we  compare  the  different  types  of  plasmonic  waveguides  that  we  have  introduced  in  this 
section,  highlighting  the  differences  among  them  and  the  potentials  that  they  may  offer  for 
different  applications. 


h.  Figure  of  merit  for  different  plasmonic  waveguides:  ID  propagation 

Following  Sections  2b-c,  ID  plasmonic  waveguides  may  be  envisioned  as  linear  arrays  of 
plasmonic  nanoparticles  in  a  transparent  background,  nanovoids  in  a  plasmonic  background,  and 
cylindrical  nanorods.  In  all  these  scenarios,  the  field  distribution  away  from  the  interface 

between  plasmonic  and  dielectric  materials  decays  as  AT  {\j/32  -  co2sbp0p j,  where  p  is  the 

radial  coordinate  with  respect  to  the  axis  of  propagation,  Kt  (.)  is  the  modified  Hankel  function 

and  1  =  1  for  azimuthally  symmetric  modes  (the  longitudinal  mode  in  the  array  geometries  and 
the  TM0  modes  in  the  cylindrical  structures),  whereas  i-  2  for  the  transverse  polarization  and 

the  hybrid  quasi- TMX  modes  [39]-[4 1  ] .  This  allows  one  to  evaluate  the  beam  cross  section  as  the 

distance  pQ  for  which  the  field  amplitude  is  decreased  to  e~l  of  its  value  at  p  =  a  .  It  is  clear  that 

p0  decreases  for  larger  Re[/l] ,  since  the  mode  is  more  concentrated  around  the  waveguide.  Not 

surprisingly,  decreasing  the  waveguide  radius  a  for  a  given  frequency,  Re  [/?]  increases,  in 

principle  without  any  limit,  supporting  the  possibility  of  highly  confined  sub-diffractive  guided 
wave  propagation. 
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On  the  other  hand,  we  may  define  the  attenuation  length  L  as  the  propagation  length  along  the 
x  axis  for  which  the  modal  field  has  decreased  to  e~'  of  its  original  value  due  to  ohmic  (and  in 

general  radiation)  losses.  It  is  noted  that,  due  to  the  same  definition  of  /? ,  L  =  ( Im  [/?]  )  ' . 


b)  Wavelength  X0  [  nm  ] 

Figure  10.  Dispersion  of  the  lateral  confinement  (a)  and  attenuation  length  (b)  for  four  different  forward-wave  1-D 
waveguides  employing  silicon-nitride  and  silver  with  a  =  25  nm  . 

As  a  first  example,  in  Fig.  10  we  compare  the  guidance  properties  of  four  different  ID  forward- 
wave  waveguides  composed  of  silver  Ag  (as  the  plasmonic  material)  and  silicon-nitride  Si3N4 
(as  the  dielectric),  in  particular  for  the  cases  of:  (a)  silver  nanoparticles  embedded  in  the  silicon- 
nitride  background  with  longitudinal  polarization,  (b)  silicon-nitride  nanovoids  in  silver 
background  with  transverse  polarization,  (c)  silver  nanorod  in  silicon-nitride  background  with 
TM0  mode,  and  (d)  a  silicon-nitride  nanocylinder  in  silver  background  with  quasi  -  TM t  hybrid 

mode.  It  is  important  to  stress  that  for  a  fair  comparison  of  the  different  performance  of  these 
different  waveguides,  in  this  section  we  have  used  the  experimental  values  of  the  silver 
permittivity  as  extracted  from  Ref.  [27].  This  is  different  from  the  simplified  Drude  model  used 
in  the  previous  section,  which  does  not  take  into  account  the  interband  transitions  of  silver  in  the 
near-UV.  We  have  preferred  the  Drude  model  in  the  previous  section  because  it  provides  a 
clearer  and  more  predictable  behavior  with  frequency,  providing  some  intuitive  insights  into  the 
effect  of  frequency  dispersion  on  the  modal  propagation  of  these  sub-diffractive  modes. 
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However,  in  this  section  we  use  realistic  experimental  values  of  silver  permittivity,  to  ensure  the 
higher  accuracy  of  these  results  even  near  the  UV  band. 


In  Fig.  10a  we  plot  the  variation  of  the  lateral  confinement  p0  versus  the  operation  wavelength 
(in  free-space)  A0 ,  whereas  in  Fig.  10b  we  report  the  corresponding  variation  of  the  attenuation 
length  L  .  In  all  these  examples  the  radius  of  the  waveguide  is  a  =  25  nm  and  for  the  array  cases 
the  distance  between  neighboring  spherical  inclusions  is  d  =  2.\ a  .  Several  features  are  evident 
from  these  plots:  first  of  all,  in  all  these  geometries  it  is  possible  to  confine  the  propagating 
beam,  over  a  given  frequency  band,  in  a  cross  section  much  smaller  than  the  usual  diffraction 
limit  (Ag/2).  These  sub-diffractive  beams  may  propagate  over  several  wavelengths  without 

diffraction  and  with  relatively  low  damping,  despite  the  strong  concentration  of  the  guided  beam 
and  the  material  absorption  (and  possibly  radiation  in  the  cases  of  transparent  backgrounds). 


Figure  1 1  -  Figure  of  merit  for  the  four  classes  of  forward- wave  ID  waveguides  in  Fig.  10  (when  available)  for  four 
different  scenarios:  (a)  silver  and  free-space  with  a  =  25  nm  ;  (b)  silver  and  silicon-nitride  with  a  =  25  nm 
(corresponding  to  Fig.  10);  (c)  silver  and  free-space  with  a  =  5  nm  ;  (d)  silver  and  silicon-nitride  with  a  =  5  nm  . 


The  available  bandwidth  is  usually  much  larger  for  cylindrical  waveguides,  whereas  the 
propagation  along  chain  of  particles  is  usually  possible  only  over  a  relatively  narrower 
bandwidth  (still  large  enough  for  many  applications)  around  the  resonance  frequency  of  the 
inclusions  composing  the  array,  which  in  this  case  happens  near  the  wavelength  for  which 
sAg  =  —2 ±leSi  N  (with  the  plus  minus  depending  on  whether  the  silver  composes  the  inclusions  or 

the  background).  Depending  on  the  frequency  band  and  the  application  of  interest,  arrays  of 
particles  or  cylindrical  waveguides  may  be  more  or  less  appealing.  For  instance,  in  this  geometry 
the  array  of  nanovoids  appears  very  suitable  in  the  optical  frequency  range,  providing  the  longest 
propagation  distance  together  with  a  good  beam  confinement.  At  lower  frequencies,  the 
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cylindrical  plasmonic  nanorod  seems  the  most  suitable  solution  for  this  combination  of  materials, 
due  also  to  the  larger  bandwidth  of  propagation. 

A  clear  general  trend  is  evident  in  Fig.  10,  which  applies  to  all  these  classes  of  plasmonic 
waveguides:  a  trade-off  should  be  sought  between  beam  lateral  confinement  and  attenuation 
length.  Even  if  in  principle  there  is  not  an  upper  limit  to  the  value  of  Re  [/?]  (and  therefore  to  the 

beam  lateral  confinement),  the  value  of  1m  [/?]  indeed  usually  grows  accordingly,  effectively 

limiting  the  possibility  of  realizing  highly  confined  ultra-sub-diffractive  waveguides  with 
reasonably  long  attenuation  lengths.  In  the  following,  therefore,  we  define  a  figure  of  merit  as  the 
ratio  F  =  L  t  pQ ,  in  order  to  effectively  compare  waveguides  with  different  sizes  and  geometries. 

Figure  11,  for  instance,  compares,  when  available,  the  four  different  scenarios  of  Fig.  10  using 
two  different  dielectric  materials  (free-space,  Fig.  lla,c;  silicon-nitride,  Fig.  llb,d)  and  two 
different  radii  (a  =  25 nm ,  Fig.  lla,b;  a=5nm,  Fig.  11c, d).  For  each  panel  only  some 
waveguides  may  support  a  propagating  mode  with  sufficient  robustness  to  propagation  despite 
the  material  absorption,  as  indicated  in  the  legend.  Even  though  the  trends  are  consistent  in  the 
four  panels,  important  differences  arise.  For  instance,  a  larger  radius  usually  allows  reaching 
larger  figures  of  merit,  due  to  the  lower  field  concentration  in  the  plasmonic  materials. 
Moreover,  the  free-space  as  dielectric  ensures  larger  figures  of  merit  at  low  frequencies,  since  the 
mode  can  be  more  spread  (with  lower  beam  concentration)  in  the  transparent  region,  even  though 
the  chain  configuration  in  the  plasmonic  background  may  achieve  relatively  large  figures  of 
merit  at  optical  wavelengths. 


Figure  12.  Similar  to  Fig.  11,  but  for  four  backward-wave  ID  waveguides  (dual  of  those  in  Fig.  11,  when  available) 
for:  (a)  silver  and  free-space  with  a  =  25  nm  ;  (b)  silver  and  silicon-nitride  with  a  =  25  nm  (corresponding  to  Fig. 
10);  (c)  silver  and  free-space  with  a  =  5  nm  ;  (d)  silver  and  silicon-nitride  with  a  =  5nm  . 
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Figure  12  refers  to  the  backward  modes  supported  by  the  dual  geometries  or  polarizations, 
employing  the  same  materials.  In  this  case,  the  modes  are  supported  only  by  some  combinations, 
since  backward  modes  are  usually  more  challenging  to  be  supported.  However,  consistent  with 
Fig.  8,  for  small  radii  the  nanovoid  linear  array  supports  a  fairly  robust  longitudinal  mode  with 
backward-wave  properties  that  may  achieve,  over  a  reasonable  bandwidth,  considerable  figures 
of  merit,  as  reported  in  Fig.  12b.  This  may  provide  some  interesting  possibility  to  realize  left- 
handed  nanotransmission  lines  at  optical  frequencies.  The  cylindrical  rods,  on  the  other  hand,  do 
not  easily  support  backward  modes  with  sufficient  robustness  to  losses.  If  the  interest  is  merely 
in  field  confinement  and  propagation  length,  the  best  option  and  the  largest  figure  of  merit  is 
usually  obtained  with  forward-mode  waveguides. 
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a)  Wavelength  X0  [  nm  ] 


Figure  13.  Similar  to  Fig.  10,  but  for  2D  geometries,  dispersion  of  the  lateral  confinement  (a)  and  attenuation  length 
(b)  for  three  different  forward- wave  2-D  waveguides  employing  silicon-nitride  and  silver  with  a  =  25  nm  . 

i.  Figure  of  Merit  for  Different  Plasmonic  Waveguides:  2D  propagation 

We  may  perform  an  analogous  analysis  for  the  different  types  of  2D  waveguides,  i.e.,  2D  planar 
arrays  of  nanoparticles  or  nanovoids  and  planar  nanolayers  with  different  combinations  of 
plasmonic  and  dielectric  materials.  In  this  scenario,  the  field  distribution  decays  in  the  transverse 
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In  the  case  of  2D  arrays  of 

inclusions,  sufficiently  small  dx  =  d  v  ensure  isotropic  propagation  in  two  dimensions,  analogous 
to  the  planar  layers.  In  the  following,  for  comparison,  we  analyze  this  situation. 

Figure  13  reports  the  dispersion  of  p0  and  L  for  three  classes  of  plasmonic  2D  forward 

waveguides,  consisting  of  a  planar  array  of  silver  nanoparticles  in  a  silicon-nitride  background 
(with  longitudinal  polarization),  a  silver  planar  nanolayer  surrounded  by  silicon-nitride  (odd 
mode)  and  a  Si3 NA  nanolayer  in  a  silver  background  (even  mode).  For  the  planar  array  of 

nanoparticles  we  have  assumed  isotropic  guidance  properties,  i.e.,  same  interspacing 
dx  -  d  =2.1  a  in  the  plane  of  propagation.  In  all  the  cases  a  =  25  nm  is  the  nanosphere  radius  or 

the  half-thickness  of  the  nanolayer.  The  dual  scenario  of  silicon-nitride  voids  in  a  silver 
background  does  not  support  a  forward  mode  with  sufficiently  low  damping  in  this 
configuration,  so  it  has  not  been  reported  in  the  figure. 


Figure  14.  Similar  to  Fig.  11,  but  for  three  forward-wave  2D  waveguides  for:  (a)  silver  and  free-space  with 
a  =  25  nm\  (b)  silver  and  silicon-nitride  with  a  =  25  nm  (corresponding  to  Fig.  8);  (c)  silver  and  free-space  with 
a  =  5  nm  ;  (d)  silver  and  silicon-nitride  with  a  =  5  nm  . 

It  is  noticed  that  the  propagation  properties  are  analogous  to  the  ID  case,  even  if  the  bandwidth 
and  guidance  properties  of  the  arrays  of  nanoinclusions  are  somehow  worsened  by  the  coupling 
between  closely  packed  parallel  arrays.  The  plasmonic  nanolayer  geometry  may  achieve  very 
long  propagation  distances,  at  the  expenses  of  worse  beam  confinement.  On  the  other  hand,  in 
the  infrared  the  dielectric  planar  gap  in  a  silver  background  performs  well,  similar  to  a  parallel- 
plate  waveguide,  confining  the  field  in  a  sub-diffractive  quasi-TEM  mode  with  long  propagation 
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distances.  In  the  2D  scenario  this  parallel-plate  silver  configuration  seems  the  most  adequate  for 
guidance  and  in  this  sense  the  advantages  of  a  plasmonic  background,  as  outlined  in  the  previous 
sections,  are  evident  in  Fig.  13. 

Figure  14  shows  the  figure  of  merit  dispersion  for  four  different  geometries,  consistent  with  the 
results  of  Fig.  11,  but  for  the  2D  case.  Here  it  is  even  more  evident  how  the  plasmonic 
background  allows  confinement  in  the  waveguide  in  a  sub-diffractive  region,  still  ensuring 
sufficiently  low-damping.  At  optical  frequencies  the  2D  arrays  of  nanovoids  may  also  be 
employed  in  this  configuration  with  good  performance.  (As  an  aside,  it  should  be  noted  that  the 
central  frequency  for  the  propagation  band  in  the  array  geometries  may  be  tailored  at  will  by 
changing  the  geometry  of  each  inclusion,  or  the  involved  materials,  as  outlined  in  the  previous 
section.  In  this  geometry  the  resonance  between  silver  and  air  or  silicon-nitride  happens  in  the 
visible,  but  THz  plasmonic  materials,  like  silicon  carbide,  or  multi-layered  particles,  may  shift 
the  propagation  band  at  will). 


Figure  15.  Variation  of  the  figure  of  merit  for  forward  ID  and  2D  waveguides  as  a  function  of  the  half-thickness  a 
for:  (a)  silver  and  free-space  at  A0  =600  nm  ;  (b)  silver  and  silicon-nitride  at  A0  =  600  nm  (corresponding  to  Fig. 
13);  (c)  silver  and  free-space  at  A0  =1.5  jum  ;  (d)  silver  and  silicon-nitride  at  A0  =1.5 /uni . 

As  a  last  example,  Fig.  15  shows  the  variation  of  F  ,  for  the  classes  of  forward  waveguides  in  ID 
and  2D  geometries,  with  the  half  thickness  a  of  the  waveguide,  for  the  same  combinations  of 
materials  as  in  the  previous  examples.  The  panels  refer  to  an  optical  wavelength  (20  =  600  nm  , 
Fig.  lOa-b)  and  a  typical  telecommunication  wavelength  ( A0  =1.5 /urn,  Fig.  lOc-d).  These  charts 
show  how  the  figure  of  merit  tends  to  increase  for  larger  waveguides  and  it  is  in  general  larger  at 
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lower  frequencies  and  for  Si3 N4  -  Ag  combinations.  The  potential  advantages  of  a  plasmonic 
background  in  guidance  and  confinement  are  evident  also  in  these  charts. 

j.  Conclusions 

We  have  considered  and  compared  in  this  section  eight  different  solutions  for  guiding  optical 
beams  with  highly  confined  sub-diffractive  cross-sections  and  reasonable  propagation  distance, 
comparing  their  potential  advantages  and  drawbacks.  In  particular,  we  have  compared,  both  in 
ID  and  2D  geometries,  plasmonic  waveguides  embedded  in  a  transparent  background,  and  the 
dielectric  waveguides  in  a  plasmonic  background,  showing  how  the  latter  may  arguably  provide 
better  figures  of  merit  in  terms  of  propagation  distances  versus  lateral  cross-section.  The 
analogies  and  differences  between  arrays  of  particles  and  continuous  rods  and  layers  have  been 
outlined,  showing  how  the  latter  may  provide  wider  guidance  bands,  whereas  the  former  may 
show  larger  figures  of  merit  and  propagation  lengths,  but  more  limited  bandwidth  of  operation 
concentrated  around  the  resonances  of  the  inclusions  composing  the  arrays.  These  results  may  be 
of  interest  for  the  design  and  realization  of  plasmonic  waveguides,  with  potentials  in  optical 
communications  and  nanocircuit  applications. 


5.  Radiation  and  Leaky-Waves  along  Linear  Arrays 

a.  Summary 

We  analyze  in  this  section  the  leaky-wave  properties  of  linear  arrays  of  plasmonic  nanoparticles. 
It  is  shown  that  such  periodic  arrays  may  support  two  orthogonal  leaky-wave  propagation 
regimes,  respectively  with  longitudinal  and  transverse  polarization.  Using  closed-form  dispersion 
relations  derived  in  the  complex  domain,  consistent  with  the  previous  sections,  we  analyze  their 
properties  in  the  leaky-wave  regime  and  we  derive  general  conditions  under  which  a  nanoparticle 
array  with  sub-wavelength  lateral  cross  section  may  support  a  radiating  leaky  mode  with 
directive  properties,  conical  radiation,  frequency  scanning  and  sufficiently  long  propagation 
distance,  paving  the  way  to  potential  applications  as  a  leaky-wave  optical  nanoantenna  with  sub- 
diffractive  properties.  Realistic  designs  and  configurations  are  presented,  considering  the 
material  dispersion  and  absorption  of  optical  materials,  for  which  we  detennine  propagation 
distance,  near-field  distribution  and  far-field  leaky-wave  radiation  pattern. 

b.  Introduction 

The  miniaturization  of  electronic  and  optical  devices  is  one  of  the  main  challenges  in  modem 
communications  and  computer  technology.  Various  concepts  and  devices,  well-established  in 
microwave  engineering,  have  been  transplanted  to  optical  frequencies,  at  which  the  characteristic 
size  and  operating  wavelength  are  orders  of  magnitude  smaller,  and  frequency  bandwidths  are 
proportionally  larger.  One  successful  example  is  represented  by  optical  nanoantennas  [74]-[82], 
which  have  been  inspired  in  recent  years  by  well-established  concepts  at  radio  frequencies[83]. 
As  another  example,  in  microwave  technology  the  electromagnetic  properties  of  periodic 
structures  play  a  crucial  role  in  several  devices.  Various  periodic  structures,  such  as  slot  arrays 
and  frequency  selective  surfaces,  are  widely  applied  as  antennas  and  filters.  Recent  advances  in 
nanotechnology  have  made  possible  to  extend  also  these  concepts  to  optical  frequencies,  where 
periodic  structures,  arrays  and  nanoscale  metamaterials  have  been  recently  investigated  for  a 
variety  of  applications. 
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As  one  of  the  interesting  applications  of  periodic  arrays  at  radio  frequencies  for  radiation 
applications,  leaky-wave  antennas  are  a  well-established  technology  that  provides  directive 
radiation  and  frequency  beam  scanning  [83]-[85].  The  recent  application  of  metamaterial 
concepts  has  provided  novel  possibilities  for  leaky-wave  antenna  design  and  operation  at 
microwave  frequencies  [86]-[88].  Translating  these  concepts  to  the  optical  domain  may  open 
new  areas  in  optical  communications,  control  of  radiation  and  optical  computing.  In  this  regard, 
periodic  arrays  of  nanoparticles  have  already  been  considered  by  various  groups  as  optical 
waveguides  with  confined  beams,  overcoming  the  optical  diffraction  limit  [89]-[98] .  Dielectric 
waveguides  are  generally  limited  by  diffraction  to  have  a  transverse  cross  section  comparable 
with  the  wavelength,  as  guided  optical  beams  tend  to  spread  in  the  background  material  when  the 
waveguide  is  too  thin  [99].  However,  the  use  of  plasmonic  materials,  and  arrays  of 
subwavelength  plasmonic  nanoparticles  in  particular  [89]-[98],  may  overcome  this  limitation  and 
confine  a  guided  optical  beam  over  a  transverse  cross-sections  significantly  smaller  than  the 
wavelength,  supporting  sub-diffractive  propagation  with  relevant  applications  in  optical 
computing  and  communications. 

This  same  nanoparticle  array,  which  is  realizable  within  available  nanofabrication  technology, 
may  also  provide  an  interesting  way  of  realizing  leaky-wave  nanoantennas  with  sub-wavelength 
lateral  cross-section  and  directive  radiation  at  a  specific  angle  in  the  far-field.  Our  group  has 
theoretically  investigated  in  the  past  guided-wave  propagation  along  linear  chains  of  plasmonic 
and  metamaterial  particles  as  optical  nanotransmission  lines  with  sub-diffractive  properties 
[100].  Our  contribution  to  this  problem  consisted  in  the  derivation  of  a  closed-form  dispersion 
relation  for  real  and  complex  dipolar  modes  supported  by  such  arrays,  with  the  only 
approximation  being  the  neglect  of  multipoles  of  higher  order  than  the  dominant  dipolar 
contribution  from  each  particle.  In  particular,  this  fonnulation  makes  it  possible  to  deal  with  the 
presence  of  realistic  losses  and  damping  for  the  guided  modes,  extending  previous  analyses  that 
were  limited  to  real  wave  numbers  to  the  complex  domain  by  an  analytic  continuation  technique 
[100].  Similarly,  this  technique  may  be  applied  to  problems  involving  radiation  losses,  coming 
into  play  when  the  wave  energy  is  not  totally  guided  along  the  particle  chain,  but  partially  leaked 
out,  as  it  happens  in  leaky  modes. 

As  mentioned  above,  the  idea  of  energy  leakage  is  widely  applied  in  microwave  engineering  to 
design  directive  radiators  with  beam  scanning  capabilities.  A  leaky  mode  is  a  fast  eigenmode  of 
the  structure  with  complex  wave  number,  whose  real  part  is  less  than  the  free-space  wave 
number  [101].  This  ensures  that  the  energy  is  not  confined  along  the  array,  and  the  Poynting 
(power  flux)  vector  points  towards  the  lateral  direction.  Provided  that  the  imaginary  part  of  the 
leaky  wave  number  is  sufficiently  small,  the  radiation  from  the  chain  may  become  very  directive, 
producing  a  conical  directive  beam  at  a  given  angle  from  the  array  axis.  At  microwaves,  leaky- 
wave  antennas  are  usually  obtained  by  perturbing  a  guided  wave  with  periodic  defects,  as  in  a 
periodically  loaded  micro-strip  line  [102]-[104].  It  is  challenging,  however,  to  produce  defects 
within  a  sub-wavelength  transverse  cross-section,  since  in  such  case  they  tend  to  weakly  interact 
with  the  mode  of  interest,  which  is  usually  weakly  confined.  This  is  another  clear  symptom  of 
the  diffraction  limit  of  guided  beams  in  free-space.  For  this  reason,  the  leaky-wave  antenna 
transverse  cross  section  is  usually  comparable  with  the  wavelength  of  operation.  In  optics, 
surface  plasmons  may  be  able  to  confine  the  energy  within  sub-wavelength  cross-sections,  and 
there  has  been  some  interest  in  using  energy  leakage  from  thin  plasmonic  films  for  near  field 
microscopy  [105]-[107]. 
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In  this  work,  stemming  from  our  previous  analysis  of  guided  modes  along  periodic  linear  arrays 
of  sub-wavelength  nanoparticles,  as  reported  in  part  in  the  previous  sections,  we  analyze  the 
potential  of  this  geometry  to  support  leaky-waves  with  directive  radiation  properties  in  the 
optical  regime,  even  in  the  limit  in  which  it  has  sub-wavelength  (i.e.,  not  limited  by  the 
diffraction  limitations  mentioned  above)  lateral  cross  section,  in  order  to  form  a  sub-diffractive 
optical  leaky-wave  nanoantenna.  This  may  lead  to  the  possibility  to  connect  distant  points  of  an 
optical  nanocircuit  board  [108]  and  create  point-to-point  links  at  the  nanoscale.  In  this  context, 
interest  in  tailoring  the  optical  radiation  from  linear  and  planar  arrays  of  nanoparticles,  forming 
Yagi-Uda  nanoantenna  arrays  [109]-[1 10]  or  planar  reflectarrays  [11 1]-[1 12],  has  been  recently 
discussed  in  several  exciting  papers.  In  the  following,  we  derive  relevant  design  parameters  and 
underline  the  fundamental  and  general  limitations  and  challenges  to  the  practical  realization  of 
leaky- wave  nanoantennas  as  linear  arrays  of  plasmonic  nanoparticles.  It  should  be  mentioned 
that  an  extensive  analysis  of  the  complex  modes  supported  by  1-D,  2-D  and  3-D  arrays  of 
magnetodielectric  particles  has  been  recently  reported  [113],  including  some  aspects  of  the 
leaky-wave  propagation  along  sub-wavelength  arrays  of  dielectric  and  magnetodielectric  spheres 
with  large  index  of  refraction.  Our  general  analysis  is  focused  here  on  plasmonic  nanoparticle 
arrays,  which  may  ensure  the  application  of  these  concepts  at  optical  frequencies  and  may 
provide  inherent  advantages  associated  with  their  anomalous  light  interaction. 


(a)  Longitudinal 
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Figure  16  -  Geometry  under  consideration:  a  linear  array  of  polarizable  nanoparticles  supporting  a  longitudinal  (a)  or 

a  transverse  (b)  eigenmode. 


c.  Theoretical  formulation 

Consider  an  infinite  linear  array  of  particles  oriented  along  the  z  axis,  periodically  located  at 
z  =  Nd ,  with  d  being  the  center-to-center  distance  and  N  being  any  positive  or  negative 
integer,  consistent  with  the  geometry  analyzed  in  Ref.  [100].  Provided  that  the  nanoparticle  size 
is  much  smaller  than  the  wavelength  of  operation,  as  assumed  in  the  previous  sections,  its  wave 
interaction  is  dominated  by  the  dipolar  scattering  and  each  element  may  be  safely  modeled  as  a 
polarizable  dipole,  fully  characterized  by  its  electric  polarizability  aee.  As  commonly  done 
[114],  and  consistent  with  the  analytical  theory  in  Ref.  [100],  if  p0  =  aee E0  is  the  dipole  moment 
induced  by  a  local  electric  field  E0  on  the  particle  at  z  =  0 ,  it  is  possible  to  derive  a  self- 
sustained  eigensolution  traveling  along  the  array  in  the  form  p  v  =  p0e  ,pNd ,  under  an  e  mt  time 
convention.  Here,  [3  is  the  complex  propagation  factor,  fully  characterizing  its  propagation  and 
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radiation  properties.  As  reported  in  several  papers  on  the  topic  [90]-[100],  the  complete 
eigenmode  spectrum  may  be  split  into  longitudinal  and  transverse  polarizations,  consistent 
respectively  with  Fig.  16(a)  and  Fig.  16(b).  The  dispersion  relations  for  these  two  polarizations 
may  be  respectively  written  [100],  consistent  with  the  previous  sections: 


L:3d 


/3  [p,  d )  -  id  f2  (p,  d )]  =  a~xe 


T :  — d 
2 


3  -r  7  r  „  /  7T  T  \  .  T 


f,  (A,  d)-idf2(p,  d)-d2ft(/],d) 


=  cc  2 


(21) 


where  d  =  k0d  ,  /?  =  /?/  k0 ,  aee  =  klaee  /  ( 6  ns{) ) ,  ktl  =  o)^£t)pit]  ,  and  s0 ,  /u0  are  the  permittivity 
and  penneability  of  background  medium,  respectively.  In  addition: 


(22) 


and  Li  v  (xj  is  the  polylogarithm  function,  as  defined  in  Ref.  [116],  Due  to  the  inherent 
periodicity  of  the  Floquet  modes  of  the  linear  chain,  we  limit  our  analysis  to  the  principal  period 

Re[/?]|  <  7r/d  . 


The  form  of  dispersion  relation  (21)  is  valid  for  any  real  or  complex  value  of  /? ,  ensuring  that  it 
may  be  employed  to  study  guided  [100]  as  well  as  leaky- wave  propagation  along  the  linear 
chains.  In  our  previous  work,  we  have  discussed  guided  propagation  along  arrays  of  extremely 
sub-wavelength  particles,  showing  that  the  condition  Im^a”1]  =  -1  is  required  for  the  involved 

nanoparticles  to  support  a  real  solution  for  ft  (guided  modes  with  no  decay)  [100].  This 
condition  is  identically  met  for  passive  dipolar  particles  only  when  absorption  may  be  neglected 
[114], [117],  as  physically  expected,  and  it  implies  that  Re[/?]>1  in  Eq.  (21).  If  the  lossless 

condition  is  not  satisfied  (Im^cf”1^]  < -1),  then  absorption  takes  place  in  the  nanoparticle  array 

and  the  eigenwave  numbers  are  necessarily  complex,  whose  imaginary  part  is  associated  with  the 
damping  caused  by  Ohmic  loss. 


Even  in  the  lossless  scenario,  however,  complex  solutions  are  allowed  when  Re[/?J<1  (fast 
leaky  modes),  when  1  <  Re <  n I d  (complex  modes)  or  when  Re[/?]|  =  Kid  (stop-band). 
In  the  following,  we  are  interested  in  leaky  modes  with  sufficiently  low  hn^/?^j,  which  may 


provide  directive  radiation  and  sustained  propagation  over  a  reasonable  electrical  length, 
analogous  to  the  operation  of  microwave  leaky-wave  antennas  [83],  For  d  <n  (sufficiently  tight 
arrays,  which  is  required  for  leaky  radiation,  as  we  note  in  the  following),  the  first-order  Bloch 
mode  dominates  the  far-field  pattern,  which  may  be  therefore  evaluated  by  simply  assuming  an 
averaged  current  line  distribution  along  the  z  axis  with  amplitude  -z<x>p0e'^z  /  d ,  consistent  with 
Ref.  [1 15],  In  this  case,  the  magnetic  potential  A  may  be  written  in  the  two  polarizations  as: 


A  = 


^AoPo 
4  d 


H, 


(i) 


(aAo2-/?2 p) 


JP- 


(23) 


35 


Approved  for  public  release;  distribution  unlimited. 


where  p  is  the  radial  coordinate  in  the  suitable  cylindrical  reference  system  with  axis  along  the 
cylinder.  The  electric  and  magnetic  far- (1  eld  distributions  may  be  easily  derived  as 


H  =  V x  A / p0,  E  =  VxH  / ) . 


This  implies  that  a  complex  value  of  /3  necessarily  requires  a  non-zero  power  flux  and  phase 
propagation  along  the  radial  direction.  In  particular,  for  Im  [/?]  sufficiently  small,  Eq.  (23) 
represents  a  standard  guided-wave  mode  for  Re[/?]>k0,  with  exponential  decay  rate  in  the 
radial  direction  given  by  Ref.  [100]: 


I:  K,()Re[/?f-*aV) 

T:  K2(^Re[/?]!-*0V) 


(24) 


and  a  leaky  mode  when  Re[/?]  <  k0 ,  with  conical  beam  radiation  at  an  angle  6  =  cos  1  |^Re[/?] 

from  the  z  axis.  In  such  case,  the  decay  rate  is  the  one  of  a  cylindrical  wave  1  /  -sfp  and  the 
corresponding  intensity  pattern  is  well  approximated  by  [84]: 

sin2  0 


m= 


(cos^-Re/?) 


tt\2 


(25) 


The  radiation  beamwidth  of  the  main  conical  lobe  is  calculated  as: 


BW  =  2  Im /?  /  sin  <90 ,  (26) 

which  ensures  that  the  directivity  of  radiation,  a  measure  of  how  oriented  and  narrow  the  far- 
field  radiation  pattern  is  towards  the  desired  direction,  is  inversely  proportional  to  hn[/?J  . 

It  should  be  emphasized  in  this  context  that  a  complex  solution  of  Eq.  (21)  implies  in  principle  a 
diverging  distribution  of  the  induced  dipole  distribution  along  the  array,  which  may  raise  some 
questions  about  the  validity  of  its  analytical  continuation  in  the  complex  domain.  It  is  noticed, 
however,  that  similar  concerns  arise  any  time  we  deal  with  complex  eigensolutions  of  the  wave 
equation,  i.e.,  in  regular  leaky- wave  configurations,  or  even  in  surface- wave  propagation  along 
lossy  interfaces.  Similar  to  such  cases,  this  divergence  arises  only  because  of  the  assumption  of 
an  infinite  array.  Leaky-wave  solutions  do  not  represent  proper  contributions  to  the  radiated 
spectrum  of  the  chain,  but  they  indeed  dominate  the  steepest-descent  approximation  in  specific 
angular  regions  of  the  visible  spectrum,  and  therefore  they  constitute  an  accurate  and  effective 
description  of  the  far-field  distribution  of  the  chain  in  a  variety  of  realistic  applications  [83].  In 
practice,  their  divergence  does  not  constitute  an  issue,  since  we  are  interested  in  solutions  with 
small  I m  [ /? ]  and  finite  chain  lengths,  for  which  the  localized  excitation  (which  may  be 

represented  by  an  emitting  molecule  or  a  quantum  dot  in  this  scenario  at  optical  frequencies)  is  at 
a  finite  location  along  the  array  [118]. 
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Figure  17  -  Variation  of  complex  /?  in  (a)  longitudinal  and  (b)  transverse  polarization  versus  the  normalized  inverse 
polarizability  of  the  nanoparticles  composing  the  array.  Here  we  consider  a  normalized  center-to-center  distance 

d  =  0.2 . 

d.  General  properties  of  the  leaky-wave  eigensolution 

In  this  section,  we  report  our  investigation  on  the  general  properties  of  the  complex  solutions  of 
Eq.  (21),  with  special  attention  to  the  leaky-wave  regime.  In  order  to  make  the  analysis  very 
general,  we  focus  in  this  section  on  the  variation  of  complex  (3  with  the  normalized  quantity 

Rc[aJ  ] ,  which  compactly  describes  the  general  properties  of  the  individual  nanoparticles 
forming  the  array.  It  is  noticed,  in  particular,  that  Im  \Jx~l  J  is  simply  associated  with  the 
absorption  properties  of  the  particles,  and  it  is  forced  to  be  Im  fa  f  J  =  -1  when  the  particles  are 
lossless.  The  available  degrees  of  freedom  to  tailor  the  leaky-wave  properties  of  the  array  are 
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therefore  compactly  represented,  to  within  the  dipolar  approximation,  by  Re[aj],  which  is  a 

function  of  the  geometrical  and  material  properties  of  the  nanoparticles.  In  the  next  sections,  we 
will  provide  specific  examples  of  nanoparticle  geometries  that  may  synthesize  the  required 
values  of  Re[a^1]  obtained  in  this  section. 


As  a  first  example,  in  Figure  17  we  report  the  variation  of  complex  /?  as  a  function  of  the 
nonnalized  parameter  d 3  Re  \jxj,  j ,  for  an  interparticle  distance  d  =  0.2  .  We  consider  here 

lossless  particles  with  Im^a”1  j  =  -1 .  As  seen  in  Fig.  17(a),  the  longitudinally  polarized 
eigenmodes  have  a  smooth  transition  from  the  guided-wave  to  the  leaky-wave  region  at 
Re[/?]  =  1.  The  lossless  nature  of  the  particles  ensures  Im  [/?]  =  0  in  the  guided  region 

Rc [/?]>•  1 .  As  the  wave  number  enters  the  region  Re [/?]<!,  the  imaginary  part  starts 


Re[/?T 


.  It  is 


increasing,  due  to  the  conical  radiation  of  the  leaky  mode  at  an  angle  Q{)  =  cos  1 

recognized  that  the  guided  modes  in  this  longitudinal  polarization  are  inherently  forward  in 
nature,  since  the  slope  3Re[/?]  /fiRe[«~1  J  is  negative.  As  explicitly  proven  in  Ref.  [100],  in 

fact,  the  slope  of  the  curves  in  Fig.  17  is  directly  related  to  whether  the  modes  are  forward 
(negative  slope)  or  backward  (positive),  which  directly  determines  the  sign  of  <3R e[/?]/<3<»  for 

passive  particles  in  regions  in  which  hn[^/?j  is  negligible. 

Also  in  the  leaky- wave  regime,  for  low  negative  slope  is  preserved,  but,  for  the  value 

aj  =  am'm  ,  the  real  part  of  ft  reaches  a  minimum  at  J3min  and  then  returns  to  Re[/?]  =  1 .  Similar 

arguments  apply  in  the  low-damping  region  aj  >  am'm  ,  ensuring  that  the  supported  longitudinal 
leaky-wave  modes  are  forward,  improper[l  19]  in  nature,  as  also  verified  by  the  fact  that 
Re[/?]/Im[/?]  >  0  ,  (phase  and  group  velocities  are  parallel  with  each  other)  in  the  region  with 

small  Im[/?J.  This  follows  from  the  modal  dependence  e,fiz ,  which  ensures  that  the  phase 
propagation  is  in  the  same  direction  as  the  power  flow  and  energy  decay,  under  the  condition 
Re  /?]/hn[/?]  >  0  . 

The  level  of  radiation  damping  monotonically  increases  with  Re^aJ  j ,  implying  that  the  range 

Rc[«J  ]  <  am jn  is  preferable  for  more  directive  radiation.  This  is  physically  expected,  since  this 
is  the  region  closer  to  the  resonance  of  the  individual  nanoparticles  composing  the  array,  which 
always  arises  at  Re^a^1]  =  0.  Longitudinal  leaky  modes  are  inherently  supported  for  positive 


values  of  Re  aj  ,  due  to  their  forward  nature,  since  causality  requires[100] 


SRcf«;'_ 


dm 


<0. 


It  is  worth  noticing  that  the  point  of  minimum  Re[/?]  =  fdmm  arises  close  to  the  crossing 
Re[/?j  =  Im[/?J  in  the  plot  of  Fig.  17a.  This  point  may  be  considered  the  cut-off  of  the  leaky- 
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wave  regime,  since  for  a~j  >  am |n  the  leaky-wave  radiation  is  damped  by  rapid  longitudinal 
decay,  and  its  directivity  is  very  limited.  The  occurrence  of  a  cut-off  for  leaky  modes  close  to 
where  Re[/?]  =  Im[/?J  is  well  known  in  a  variety  of  leaky- wave  antennas[85],  and  it  is  verified 

in  this  geometry  for  different  values  of  d  in  Fig.  18.  It  is  interesting  to  note  that  this  cut-off 
arises  here  around  the  region  of  minimum  ft  . 

In  the  transverse  polarization  (Fig.  17b),  the  guided  confined  branch  (right  in  the  figure)  is 
inherently  backward  in  nature,  since  5Re[/?]/fiRe[ag~1]  >  0 ,  due  to  similar  arguments  as 

outlined  above,  and  consistent  with  analogous  findings  in  thin  plasmonic  films  [  120]-[  121]  and 
optical  metamaterials  [122].  As  outlined  in  Ref.  [100],  a  second,  less  confined,  forward  branch  is 
also  present  in  the  guided  regime,  of  less  interest  from  the  practical  point  of  view,  since  it  is  very 
similar  to  a  plane  wave  traveling  unperturbed  in  the  background  with  very  limited  confinement. 
A  complex  branch  stems  from  the  contact  point  between  these  two  guided  modes,  which  enters 
the  leaky-wave  regime  for  sufficiently  negative  Re  [a”1!.  The  dispersion  of  Re[/?J  with 

frequency  in  this  regime  decreases  monotonically  from  +1  to  -1,  for  decreasing  Re^erJ  ] , 

crossing  the  axis  Re[/?]  =  0.  For  this  specific  value  of  inverse  polarizability,  the  leaky  mode 

passes  from  backward  proper  [123]  (for  less  negative  Re  [a”1])  to  forward  improper  operation. 
It  is  evident  that  in  this  polarization  we  are  mostly  interested  in  the  backward  region,  which 
ensures  smaller  damping  factor  Im  J .  As  expected,  also  in  this  polarization  the  most 
interesting  region  arises  closer  to  the  resonance  of  the  individual  nanoparticles,  i.e.,  here  for  less 
negative  values  of  Re[aJ  ] .  The  leaky- wave  branch  is  connected  to  the  guided  branches 

through  a  complex  modal  regime,  which  is  typical  of  a  transition  between  leaky-wave  modes  and 
a  two-branch  guided- wave  regime  [85].  In  this  transition  region,  the  mode  does  not  radiate  and 
propagates  with  complex  wave  number,  whose  real  part  is  very  close  to  the  one  of  free-space, 
and  non-zero  imaginary  part. 

It  is  interesting  to  stress  that  the  inherent  backward  propagation  of  guided  and  leaky-wave  modes 
with  transverse  polarization  may  be  appealing  in  the  framework  of  negative-index  propagation, 
and  this  guided  regime  has  been  exploited  to  realize  double -negative  metamaterials  in  the  visible 
[120]-[122],  In  terms  of  leaky-wave  radiation,  backward  radiation  may  be  of  interest  to  increase 
the  degrees  of  freedom  in  tailoring  and  directing  the  optical  radiation,  but,  as  we  show  in  the 
following,  it  is  intrinsically  less  efficient  than  the  forward  longitudinal  mode.  Farther  from 
resonance,  outside  the  leaky-wave  regime,  both  polarizations  have  a  stop-band  region  with 

Re[yS]=l,  in  which  the  imaginary  part  grows  in  magnitude,  the  propagation  is  evanescent  in 

nature  and  the  damping  is  significantly  large.  In  the  following,  we  analyze  more  in  detail  the 
dispersion  of  the  leaky-wave  modes  as  a  function  of  the  interparticle  distance  d  and  of  the 
nanoparticle  polarizability,  with  the  goal  of  optimizing  the  leaky-wave  radiation  in  the  two 
polarizations,  and  of  analyzing  the  fundamental  limitations  and  possibilities  of  leaky  radiation  at 
optical  frequencies. 
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i.  Longitudinally  polarized  modes 

Consistent  with  Fig  17(a),  Re[/?]  =  1  constitutes  the  boundary  between  guided- wave  and  leaky- 

wave  operation  for  longitudinal  polarization.  Formally,  the  leaky- wave  regime  is  bounded  by  the 
following  conditions  on  the  nanoparticle  inverse  polarizability: 


%{?>)  + Cl,(2d)  +  dCl2(2d) 


d 


< 


lRe  [a«] 


<  Re 


fi{P>d)-idf2[fi,d) 


/SeC_  - 
Re[^]=l 


(27) 


where  fN[/3,d^  are  defined  in  (22),  ClN  (0)  are  the  Clausen’s  functions  [116],  which  are  real 
for  real  argument  and  Q (.)  is  the  Riemann  Zeta  function.  The  left-hand  side  has  been  written  in 
closed-form  using  the  properties  of  the  polylogarithm  functions  for  real  argument: 


Lil(ew)  =  Cll(0)  +  i 


( n~° ) 


Li2(ew)  =  - — e^ln  e\iCl1{6)  O<0<2tt 


Li3(ew)  =  CI3(0)  +  i 


.  0{n -6){2n -6 ) 


12 


(28) 


Figure  18  -  Guided-  and  leaky-wave  regions  for  longitudinal  polarization.  The  solid  blue  and  dashed  red  curves  are 
respectively  the  loci  of  real  solutions  ft  =  n /d  and  P  =  1 ,  which  delimit  the  guided-wave  regime.  The  dotted  green 

line  defines  the  upper  limit  of  the  leaky-wave  regime.  The  black  dots  denote  the  locus  Re|  /t]  =  /?mln  ,  which  may  be 

considered  the  cut-off  of  the  leaky-wave  regime. 
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Figure  18  shows  a  map  of  the  different  ranges  of  guided- wave,  leaky-wave  radiation,  or  stop- 
band,  as  a  function  of  d  and  Re  [a"1].  The  dashed  red  line  in  this  plot  represents  the  locus 

(3  =  1,  which  separates  the  guided-wave  propagation  (below)  and  leaky-wave  radiation  (above). 
The  dotted  green  line  represents  the  upper  boundary  of  the  leaky-wave  regime,  for  which  /?  e  C 
and  Re[/?]  =  1.  In  this  plot,  we  have  also  considered  the  locus  Re[/?min]  (black  dotted  line), 

which  may  be  considered  the  cut-off  for  leaky-wave  propagation,  as  discussed  above.  The  solid 
blue  line  corresponds  to  /?  =  n  /  d ,  which  is  the  lower  boundary  of  the  guided  regime.  The 
regions  above  the  leaky-wave  region  and  below  the  guided-wave  region  are  stop-band  regions, 
where  modes  decay  very  fast  along  z  ,  and  are  not  of  interest  for  guidance  or  radiation  purposes. 


Re[p] 

Figure  19  -  Variation  of  the  ratio  £  =  Re[/?]/lm[/?J  for  the  supported  leaky-wave  modes  of  the  nanoparticle  chain 
of  Fig.  16  for  longitudinal  polarization,  varying  the  center-to-center  distance. 

It  is  seen  how  all  the  boundary  curves  converge  at  d  =  n  which  represents  the  maximum 
interparticle  distance  for  supporting  guided  or  leaky  modes  along  arrays  of  sub-wavelength 
nanoparticles.  Moreover,  the  leaky-wave  region  widens  up  around  d  =  2  and  it  is  centered  above 
the  resonance  condition  for  the  individual  nanoparticles  Rej^a^,1  J  =  0 .  In  the  limit  <7  — >  0  ,  the 

leaky- wave  range  Eq.  (27)  tends  to  a  single  point  with  value  <7 3  Re[at~1]  =  6^(3)  -  7.21 , 
implying  that  too  closely  packed  chains  provide  a  very  limited  leaky-wave  radiation  bandwidth. 

One  of  the  relevant  figures  of  merit  for  leaky  modes  is  the  ratio  %  =  Re[/?]/hn[/?J .  A  lower 
Re[/?]  may  be  desirable  for  radiation  closer  to  the  normal  to  the  array,  but  this  is  usually 
accompanied  by  a  larger  hn[^/?J,  which  implies  shorter  propagation  distance,  and  inherently 
lower  directivity.  As  mentioned  above,  the  cut-off  of  the  leaky  mode  may  be  defined  by  £  =  1 . 
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Overall,  a  larger  value  of  E,  ensures  larger  directivity  and  radiation  farther  from  the  array  axis, 
both  desirable  features  of  a  leaky- wave  antenna. 

Figure  19  shows  the  variation  of  log  c  versus  Re[/?J  for  different  values  of  interparticle 
distance.  The  ratio  E,  tends  to  infinity  for  Re[/?]-l,  since  we  are  operating  near  the  guided- 

wave  regime  and  lossless  particles  are  being  considered  here.  This  region  is  characterized  by 
end  lire  radiation,  consistent  with  the  limit  of  a  surface  mode  propagating  along  the  chain.  A 
wider  range  of  Re[/?]  implies  that  energy  may  be  coupled  into  a  broader  angular  spectrum, 
which  is  more  appealing  for  antenna  applications.  Fig.  19  confirms  that  better  ratios  £  and  wider 
variation  along  Re[/?]  may  be  obtained  by  choosing  a  smaller  value  of  d  .  This  is  to  be 

expected,  since  the  nanoparticles  in  this  regime  are  tightly  coupled,  ensuring  more  flexibility  in 
tenns  of  guidance  and  radiation.  Consistent  with  Fig.  17,  however,  the  available  bandwidth  of 
leaky-wave  operation  shrinks  down  for  smaller  values  of  d  .  It  should  be  stressed,  in  addition, 
that  small  interparticle  distance  necessarily  requires  nanoparticles  with  smaller  diameters,  which 
in  turn  implies  higher  individual  Q  factors  and  more  sensitivity  to  losses.  We  discuss  these 
aspects  in  the  following  section,  when  we  consider  specific  models  for  the  nanoparticle 
geometry. 


ii.  Transversely  polarized  modes 

Transversely  polarized  leaky  modes  behave  quite  differently.  As  discussed  above,  guided-wave 
and  leaky-wave  regions  are  separated  by  a  complex  transition  region,  not  present  in  the 
longitudinal  polarization.  By  setting  Re[/?]  =  ±1  and  solving  for  the  corresponding  hn[/?]  in 

Eq.  (21),  we  can  obtain  the  range  of  polarizability  values  that  support  leaky-wave  propagation  in 
this  regime.  This  condition  may  be  formally  expressed  as: 


>3  (A d)-id  f2(p,d)-d2 f\p, d) 


Re[/J]=-1 


fi{p,d)-id  f2(p,d)-d2fx(/3,d) 


Re[/?]=1 


(29) 


where  fN[fi,d^  are  defined  in  (22). 


Figure  20  shows  the  different  modal  regions  for  transversely  polarized  modes  as  a  function  of  d 
and  a~l,  analogous  to  Fig.  18.  Like  the  longitudinal  case,  the  leaky-wave  regime  converges  to  a 
single  point  for  <7  — >  0  ,  implying  that  also  in  this  polarization  the  leaky-wave  regime  is  of 
narrow  bandwidth  for  very  tight  nanoparticles.  On  the  other  extreme,  towards  d  =  n  ,  the  modal 
region  widens,  ensuring  a  relatively  broad  range  of  normalized  polarizability  values  that  support 
leaky  modes.  As  we  will  point  out  in  the  following,  however,  the  corresponding  Im[/?]  is  rather 

large  for  this  range  of  interparticle  distance. 

Figure  21  shows  the  variation  of  E,  with  the  array  properties  in  this  polarization,  analogous  to 
Fig.  19.  It  is  evident  comparing  the  two  figures  that  it  is  more  challenging  to  obtain  a  reasonably 
large  figure  of  merit  in  this  polarization.  As  anticipated  earlier,  the  region  of  most  interest  is 
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localized  in  the  backward-wave  region,  for  positive  Re[/?]  (right  portion  of  Fig.  21).  Due  to  the 

presence  of  a  complex  transition  region  between  guided-wave  and  leaky-wave  modes,  in  this 
polarization  the  imaginary  part  is  not  negligible  even  for  values  Re[/?]  - 1“,  and  the 

ratio  £  is  never  remarkably  large.  These  results  confirm  the  general  dispersion  properties  of 
transverse  modes  highlighted  in  Fig.  17b. 


Figure  20  -  Analogous  to  Fig.  18,  guided-  and  leaky-wave  regions  for  transversely  polarized  modes. 
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Figure  21.  (Color  online)  Analogous  to  Fig.  19,  variation  of  E,  vs.  Re[/?]  for  transversely  polarized  leaky  modes, 

varying  the  center-to-center  distance. 
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Before  concluding  this  section,  it  should  be  mentioned  that  some  of  the  general  modal  features 
highlighted  here  for  linear  arrays  of  nanoparticles  may  be  obtained  in  thin  plasmonic  films,  in  the 
limit  in  which  the  array  density  increases.  This  is  consistent  with  recent  analyses  of  complex 
modes  along  such  geometries  [106],  which  have  also  highlighted  the  presence  of  backward-wave 
propagation  [120]-[121]  for  transverse  polarization,  consistent  with  the  general  results  presented 
here. 


Figure  22  -  (a)  Guided-wave  and  leaky-wave  regions  for  longitudinal  polarization,  as  a  function  of  nanosphere 
permittivity  and  interparticle  distance.  The  guided-wave  regime  is  supported  between  the  bold  lines,  while  the  leaky- 
wave  region  is  bounded  by  thinner  lines,  (b)  Loci  of  constant  £  =  Re[/?]/lm[/?  J  in  the  leaky-wave  region  for 
r)  =  2.2  .  Blue  solid  lines  delimit  the  guided-wave  and  leaky-wave  regions. 
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e.  Realistic  models  for  nanoparticles 

In  the  previous  section,  we  have  analyzed  the  general  conditions  and  limitations  for  leaky-wave 
propagation  along  sub-wavelength  nanoparticle  chains,  considering  a  general  model  for  the 
nanoparticle  polarizability.  In  particular,  we  have  shown  that  longitudinal  forward  leaky-wave 
modes  may  provide  better  directivity  properties  than  transverse  modes,  due  to  their  significantly 
larger  value  of  £  for  the  same  interparticle  distance.  Moreover,  we  have  outlined  the  range  of 

Re  [a”1]  required  to  sustain  leaky- wave  radiation.  In  this  section,  we  will  consider  realistic 

nanoparticle  geometries  to  apply  the  previous  general  results  to  several  practical  designs  for 
optical  leaky-wave  arrays. 

In  practice,  Re  [a"1]  is  determined  by  the  specific  design  and  shape  of  the  nanoparticles  forming 

the  array.  In  the  case  of  spherical  nanoparticles  of  radius  a  and  pennittivity  s  ,  for  instance,  we 
obtain  [100]: 

Re[a.;']  =  |(^)-5^'^,  (30) 

L  J  2  £-£0 

in  the  quasi-static  limit  of  interest  here.  Other  possible  geometries  of  interest  may  be  represented 
by  coated  spheres,  with  permittivities  £x  and  s2  and  ratio  of  inner  to  outer  radius  y  ,  for  which: 


Re[a«]  =  |(V)"3 


2 Y  if 2  £’i)(£-2  ^0)  (^i  ^^2)  (^2 

Y  {£2~£l)  (2^2  +  £0  )  —  (^l  +  2S2  )(^’2  —  £0  ) 


(31) 


or  nanodiscs  of  radius  a  ,  height  y  a  and  permittivity  s ,  for  which  the  transverse  polarizability 
is  [124]: 


1  1  (fzZ. 

Mi-rfl  r 


-  arctan 


r 


y 


2/ y 

£  /  £0  - 1 


(32) 


It  is  evident  that  there  is  a  wide  range  of  flexibility  in  the  shape  and  material  properties  of  the 
nanoparticles  to  tailor  the  value  of  Re  [a”1]  at  the  frequencies  of  interest.  In  the  following,  we 

focus  on  homogeneous  nanospheres  [Eq.  (30)]  and  analyze  how  their  design  parameters  affect 
the  leaky-wave  dispersion.  Analogous  results  may  be  derived  for  different  shapes  and 
geometries. 

Before  analyzing  in  detail  the  nanosphere  problem,  it  is  relevant  to  highlight  a  common  trend  in 
the  previous  formulas  (30)-(32):  as  expected,  the  value  of  Re[ag“1]  tends  to  diverge  for  small 

nanoparticles  (k0a)  — >  0  .  This  is  to  be  expected,  since  a  small  nanoparticle  is  usually  very  far 

from  its  individual  resonance  Re  \jxf  J  =  0 .  On  the  other  hand,  leaky-wave  radiation  requires 

finite  values  of  Re ] ,  as  shown  in  the  previous  section.  This  implies  that  the  operation  of 

these  leaky-wave  nanoantennas  with  sub-diffractive  lateral  cross  section  will  arise  in  the 
frequency  range  for  which  the  numerator  in  the  right-hand  side  of  Eqs.  (30)-(32)  is  close  to  zero, 
i.e.,  near  a  plasmonic  resonance  for  the  specific  shape  of  interest.  For  larger  d  this  condition 
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becomes  more  and  more  stringent,  since  Re[or(J  J  is  required  to  be  closer  to  zero.  This  is 

reflected  in  a  general  trade-off  between  size  of  these  leaky-wave  antennas  and  their  bandwidth 
and  robustness  to  the  presence  of  loss  and  disorder.  We  discuss  these  aspects  in  further  detail, 
specifically  applied  to  spherical  nanoparticles,  in  the  following. 

i.  Leaky-wave  modal  dispersion  with  the  nanosphere  permittivity 

For  spherical  nanoparticles,  we  may  use  Eq.  (30)  to  determine  the  range  of  permittivities  s  that 
may  allow  leaky- wave  propagation  along  the  nanoparticle  chain.  Figure  22a  shows  the 
longitudinal  guided  and  leaky  modal  regions  in  the  diagram  of  s  !  £0  vs.  d  ,  for  different  values 
of  the  nanosphere  radius  a  .  The  different  curves  refer  to  different  ratios  ij  =  d  /  a  and  we  have 
used  shadowing  to  highlight  the  guided-wave  and  leaky-wave  regions  in  the  case  7  =  3  .  As  it 
may  be  seen,  the  leaky-wave  region  requires  more  negative  pennittivity  values  than  the  guided- 
wave  region,  which  is  centered  at  the  resonance  condition  of  the  individual  nanospheres 
s  =  -2 s0 .  Denser  chains  support  a  wider  range  of  permittivities  to  achieve  leaky-wave 
propagation,  since  the  pennittivity  range  gets  wider  for  smaller  values  of  77  (of  course  there  is  a 
geometrical  limit  of  r/  >  2  to  consider  in  the  design).  This  is  reflected  in  wider  bandwidths,  as 
negative  permittivity  is  necessarily  dispersive  with  frequency  [125].  Consistent  with  the  previous 
section,  in  the  mathematical  limit  d  — »  0  leaky  waves  are  not  supported,  but  the  permittivity 
region  rapidly  widens  up  for  slightly  larger  values  of  d  .  Fig.  22b  shows  the  loci  of  constant 
£  =  Re[/?]/lm[/?]  for  rj  =  2.2,  as  an  example.  Consistent  with  the  results  of  the  previous 

section,  it  is  seen  that  low  attenuation  rate  is  achieved  close  to  the  boundary  of  the  guided-wave 
mode  region,  corresponding  to  end-fire  radiation.  However,  relatively  large  values  of  £  may  be 
achieved  even  farther  away  from  the  guided-wave  regime,  which  may  provide  conical  radiation 
off-axis.  Moreover,  the  natural  pennittivity  dispersion  of  metals  may  provide  frequency  scanning 
for  the  conical  beams  radiated  by  the  chain. 

Figure  23  shows  analogous  plots  for  the  transverse  polarization.  Due  to  the  backward  nature  of 
guided  and  leaky  modes  in  this  polarization,  less  negative  values  of  pennittivity  are  required  as 
compared  to  the  guided  regime.  Also  in  this  case,  by  decreasing  the  value  of  77  the  leaky-wave 
operation  broadens  in  bandwidth.  Comparing  Figs.  22  and  23,  it  is  seen  that  longitudinal  leaky 
modes  have  a  broader  leaky-wave  region  and  comparatively  larger  values  of  E,  ,  implying  that 
they  may  outperform  the  transverse  backward-wave  leaky  modes  in  terms  of  directivity  and 
bandwidth  of  operation.  These  results  are  consistent  with  the  discussion  in  the  previous  section, 
but  applied  here  specifically  to  the  nanosphere  geometry. 

ii.  Effects  of  absorption  and  material  loss 

In  this  section  we  relax  the  assumption  that  material  absorption  and  losses  are  negligible  in  the 
materials  composing  the  chain,  i.e.,  Im[a“1j<-1.  This  is  a  relevant  aspect  to  consider,  since 

negative  pennittivity,  required  to  support  subdiffractive  leaky-wave  operation,  is  usually 
combined  with  finite  absorption  [125].  Material  losses  are  known  to  play  a  relevant  role  in 
plasmonic  devices  with  subwavelength  cross  sections,  such  as  nanoparticle  waveguides  [100], 
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[126]-[127].  In  the  case  of  lossy  materials,  the  quasi-static  inverse  polarizability  is  related  to  the 
complex  permittivity  s  =  sr  +  iat  as: 


Re  [a«] 
lm[C] 


_3  / ,  \-3  (g;.  +2g0)(gr  g0)  +  g,- 

dV)  /  N.2  2 

2  +  £) 


-1 - (hflf- 

j  (sr-soy+sf 


(33) 


For  low-loss  particles,  of  interest  here,  st  is  small  and  the  associated  additional  contribution  to 
Im[aee1]  provides  a  first-order  perturbation  of  the  lossless  results  derived  above. 


d 


d 


Figure  23  -  Analogous  to  Fig.  22,  but  for  transverse  polarization. 
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Figure  24  shows  the  dispersion  of  Re[/?]  and  Im[/?J  versus  sr  for  longitudinally  polarized 

modes,  for  d  =  0.1,  ?j  =  d  /  a  =  2.\  and  different  levels  of  material  absorption  <y  .  It  is  interesting 

to  see  how  in  the  guided- wave  region  a  moderate  increase  of  st  principally  affects  1m  Jd  J ,  as 

expected,  but  leaves  unaltered  Re[/?]  and  correspondingly  the  phase  velocity.  Since  the 

transition  towards  the  leaky-wave  regime  is  continuous  for  this  polarization,  the  presence  of 
material  loss  implies  a  reduction  of  the  achievable  values  of  % ,  even  near  the  guided-wave 

region.  In  the  leaky- wave  region,  however,  the  trend  is  opposite:  Im[/?J  is  not  sensibly  altered, 
being  mainly  dominated  by  radiation  losses  (the  mode  is  less  confined  to  the  particles),  and  the 
additional  small  loss  mainly  affects  the  angle  of  radiation  and  Re  [/?  J  .  It  is  worth  noticing  that  a 

complex -valued  transition  region  may  arise  for  relatively  larger  values  of  <y ,  for  which 

Re  |>]  =  1. 


-7.0  -6.8  -6.6  -6.4  -6.2  -6.0  -5.8  -5.6  -5.4  -5.2 

S 

r 

Figure  24  -  Variation  of  Re \p"\  and  Im[/?J  for  <7=0.1  and  rj  =  d  /  a  =  2.\  in  the  longitudinal  polarization 

regime,  varying  the  imaginary  part  of  permittivity.  The  inset  plot  shows  a  zoom  in  the  transition  region  for  the  case 

£,  =0.1. 

A  zoom  of  this  transition  region  for  si  =0.1  is  reported  in  the  inset  of  Fig.  24.  The  figure 

confirms  that  realistic  levels  of  absorption  in  optical  materials  may  provide  the  possibility  to 
realize  nanoantennas  with  sub-diffractive  lateral  cross  section  able  to  sustain  such  longitudinal 
leaky  modes  with  directive  radiation  properties. 

Figure  25  shows  the  variation  of  £  for  increased  material  absorption  in  the  case  of 
rj  =  d  /  a  =  2.1  and  d  =  0.1 ,  both  for  longitudinal  (a)  and  transverse  (b)  polarization.  The  loss 
effect  is  more  evident  in  the  longitudinal  case,  since  the  transverse  polarization  has  much  larger 
radiation  losses.  Still,  the  levels  of  and  correspondingly  of  directivity,  achieved  in  the 
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longitudinal  polarization  remain  substantially  larger  than  in  the  transverse  case,  even  after 
considering  realistic  absorption  levels. 


Re[P] 


Figure  25  -  Variation  of  c  for  different  levels  of  material  loss  c  in  longitudinal  (a)  and  transverse  (b)  polarization. 
In  both  plots,  we  have  considered  f)  =  d  t  a  =  2.1  and  d  =0.1. 

Hi.  Realistic  plasmonic  materials 

The  results  of  the  previous  subsection  imply  that  chains  of  metamaterial  or  plasmonic 
nanoparticles  with  negative  permittivity  and  moderate  losses  may  provide  a  promising  mean  to 
realize  a  leaky- wave  nanoantenna  with  sub  wavelength  transverse  cross-section.  For  this  purpose, 
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noble  metals,  or  combinations  of  noble  metals  and  dielectrics,  may  be  chosen  to  realize  such 
nanochains,  following  the  design  guidelines  represented  by  Eqs.  (30)-(32). 


Figure  26  -  /3d  vs  k0d  diagrams  and  c  for  longitudinal  modes  supported  by  silver  arrays  with  i)  =  2.1. 
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Metallic  nanoparticles  made  of  silver  or  gold,  for  instance,  have  shown  moderate  guidance 
properties  in  the  optical  regime  [93].  In  this  subsection,  we  consider  the  realistic  properties  of 
noble  metals  in  the  realization  of  these  nanoantennas.  For  simplicity,  we  focus  on  nanospheres 
and  on  longitudinally  polarized  modes,  which  ensure  better  radiation  performance  and  more 
robustness  to  the  effect  of  absorption  in  the  materials  under  consideration. 


Figure  27  -  fid  vs  k0d  diagrams  for  longitudinal  modes  supported  by  dielectric  arrays  with  i)  =  2.1 . 
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Figure  26  shows  the  complex  dispersion  relations  /3d  -  k0d  for  linear  arrays  composed  of  silver 

nanospheres,  considering  experimental  values  of  permittivity  [128],  frequency  dispersion  and 
loss.  In  this  case,  we  have  chosen  77  =  2.1  and  nanosphere  radii  of  25  (blue  solid  line)  and  50  mn 
(red  dashed).  Fig  26a  reports  the  real  part  of  J3d  and  Fig.  26b  the  corresponding  imaginary  part. 
It  is  seen  that  larger  particles  may  provide  wider  bandwidth  of  leaky-wave  radiation,  due  to  the 
inherently  larger  period,  and  they  are  inherently  more  robust  to  the  presence  of  loss,  consistent 
with  analogous  results  in  the  guided  region  [100].  The  shadowed  regions  in  the  figure  indicate 
these  guidance  regions. 

Fig.  26c  reports  the  calculated  values  of  £  for  the  leaky-wave  operation.  Significantly  large 
values  may  be  achieved  near  the  endfire  radiation,  despite  the  presence  of  loss  and  the  overall 
sub-diffractive  lateral  cross-section  of  these  nanoantennas.  These  results  are  particularly 
encouraging  for  the  realization  of  these  concepts  using  arrays  of  subwavelength  silver 
nanoparticles. 


iv.  Comparison  with  dielectric  nanosphere  arrays 

The  previous  results  imply  that  plasmonic  nanoparticles  may  represent  a  promising  means  for  the 
realization  of  sub-diffractive  leaky- wave  nanoantennas.  In  this  subsection  we  compare  their 
performance  with  the  one  of  dielectric  nanoparticles,  focusing  in  the  range  d  <n .  Complex- 
wave  propagation  along  arrays  of  spheres  with  large  values  of  constitutive  parameters  has  been 
considered  in  Ref.  [113].  Figure  27  shows  the  dispersion  of  complex  modes  along  a  dense  array 
( 77  =  2. 1 )  of  dielectric  spheres  with  s  =  5^0  and  with  e  =  45^0 .  For  consistency  with  the  previous 

results,  we  show  only  guided  modes  supported  by  the  induced  electric  dipoles  along  the  array  for 
longitudinal  polarization,  although  for  large  dielectric  constants  magnetic  modes  are  also 
available.  Fig.  27a  and  27b  report  the  dispersion  diagrams  for  Re[/?]  and  the  corresponding 

Im  [/?] ,  respectively.  For  the  low  pennittivity  spheres  (blue  line),  guided  modes  are  not 

available  in  this  low-frequency  regime,  as  expected,  and  a  small  complex  branch  is  visible  near 
the  light  line.  Since  we  are  far  from  resonance,  however,  the  value  of  £  is  always  less  than  unity, 
implying  poor  radiation  properties,  as  expected.  Drastically  increasing  the  nanosphere 
pennittivity  it  is  possible  to  induce  electric  dipole  resonances,  despite  the  sub-wavelength  size  of 
the  particles.  In  this  situation,  guided-wave  and  leaky-wave  regimes  are  available,  and  the 
dispersion  diagrams  are  characterized  by  nanow  guided-wave  regions  (highlighted  by  the 
shadowed  regions)  connected  by  leaky-wave  branches.  In  some  frequency  ranges,  significant 
directivity  may  be  achieved,  although  the  radiation  is  limited  to  grazing  angles,  close  to  the  light 
line  in  the  diagrams  of  Fig.  27.  It  is  evident  that  large  permittivity  spheres  may  be  also  effective 
in  supporting  sub-diffractive  leaky  radiation,  although  the  efficiency  and  directivity  values 
achieved  in  this  example  are  lower  than  for  plasmonic  particles  and  it  may  be  challenging  to 
realize  such  large  values  of  permittivity  at  visible  wavelengths.  Plasmonic  materials  with  the 
required  values  of  permittivity,  on  the  contrary,  are  naturally  available  at  these  frequencies,  and 
their  dispersion  may  naturally  provide  a  larger  degree  of  frequency  scanning  compared  to  large 
permittivity  materials.  We  discuss  these  features  further  in  the  next  section. 
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Figure  28  -  Magnetic  field  and  power  flow  distribution  for  a  nanoparticle  chain  operating  in  the  leaky-wave  regime 
[(a)  and  (c),  at  690nm  wavelength]  and  in  the  guided  propagation  regime  [(b)  and  (d),  at  600nm]. 


f.  Full-wave  numerical  simulations 

In  the  previous  sections,  we  have  used  the  analytical  fonnulation  (21)  to  derive  the  fundamental 
properties  of  leaky-wave  propagation  and  radiation  along  infinite  arrays  of  subwavelength 
plasmonic  nanoparticles  to  within  a  dipolar  approximation.  In  this  section,  we  validate  the 
previous  analytical  model  by  simulating  realistic  finite  arrays  of  silver  nanoparticles  with  finite- 
integration  technique  commercial  software  [129],  in  order  to  determine  the  radiation  patterns  of 
such  leaky  modes  in  a  practical  realization,  considering  also  the  complete  multipolar  coupling 
among  closely  spaced  nanoparticles. 
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In  our  numerical  simulations,  we  have  fixed  the  particle  size  to  a  —  50nm ,  center-to-center 


distance  d  =  110nm,  and  we  have  used  a  Drude  permittivity  model  s  =  8x 


f(f-ir) 


,  with 


ex  =  5.0 ,  /  =  2175  THz  and  y  =  4.35  THz  ,  which  describes  with  good  approximation  the  silver 

dispersion  in  the  range  of  frequencies  of  interest  [130].  The  overall  length  of  the  chain  is 
L  -  7  pm,  sufficiently  long  to  ensure  that  significant  part  of  the  power  coupled  to  the  leaky  mode 
has  been  radiated.  The  array  is  excited  by  an  optical  source  (i.e.,  an  emitting  molecule  or  a 
quantum  dot)  longitudinally  polarized  along  the  array  axis,  to  ensure  proper  coupling  with  the 
longitudinal  leaky  modes  supported  by  the  array. 


We  have  verified  in  our  simulations  that  the  dispersion  of  leaky-wave  and  guided  modes  along 
the  array  is  qualitatively  consistent  with  our  analytical  predictions.  Clearly,  the  nature  of  our 
analytical  technique  neglects  higher-order  multipolar  coupling  between  the  closely  spaced 
nanoparticles,  which  is  reflected  in  a  quantitative  difference  in  the  prediction  of  the  frequency 
range  for  leaky-wave  radiation,  but  qualitatively  the  results  are  in  good  agreement  with  the 
previous  sections.  As  an  example,  Figure  28a  reports  the  nonnal  magnetic  field  distribution  at 
the  operating  wavelength  A0  =  690  nm,  which  is  in  the  leaky -wave  regime  for  this  array. 

Similarly,  Fig.  28b  reports  the  corresponding  distribution  at  A0  =600nm,  for  which  the  chain  is 

in  its  guided  regime.  It  is  evident  that  the  pennittivity  dispersion  of  silver  allows  tuning  the 
guidance  properties  of  the  supported  mode  from  a  slow  mode  with  short  guided  wavelength,  as 
in  Fig.  28b,  confined  along  the  structure,  to  a  much  faster  mode,  which  produces  leaky-wave 
radiation  in  free-space  with  conical  directive  properties.  The  difference  in  phase  velocity 
between  the  two  simulations  is  striking,  considering  that  the  free-space  wavelength  difference 
between  the  two  cases  is  only  15%,  and  it  is  consistent  with  our  analytical  theory.  Away  from 
the  chain,  the  leaky-wave  (Fig  28a)  couples  to  free-space  radiation,  drastically  different  from  the 
guided  propagation  in  Fig.  28b,  which  decays  exponentially  far  away  from  the  chain  axis.  The 
leaky-wave  far-field  extends  laterally  and  propagates  with  oblique  wave  fronts,  consistent  with 
the  previous  analytical  results.  Figs.  28c  and  28d  show  a  zoom  in  the  dashed  regions  of  the  two 
panels  of  Figs.  28a  and  28b,  reporting  the  power  flow  (Poynting  vector)  distribution.  The  power 
flow  shows  significant  lateral  energy  leakage  in  the  leaky-wave  scenario  of  Fig.  28c.  In  contrast, 
at  the  wavelength  A()  =600nm  (Fig.  28d),  the  power  flow  is  confined  and  bounded  parallel  to 

the  chain,  rapidly  decaying  away  from  its  axis.  It  is  remarkable  that  these  full- wave  results 
qualitatively  confirm  with  very  good  precision  the  analytical  results  in  the  previous  sections,  and 
in  particular  the  possibility  to  create  a  leaky-wave  nanoantenna  composed  of  subwavelength 
nanoparticles  composed  of  realistic  plasmonic  materials.  The  fact  that  our  full-wave  results  take 
into  account  the  full  coupling  among  the  neighboring  particles,  and  not  just  their  dipolar 
(dominant)  contribution,  slightly  shifts  the  guidance  and  leaky-wave  frequency  ranges  from  our 
analytical  predictions  in  Fig.  26,  but  qualitatively  these  results  confirm  the  possibilities  noted  in 
the  previous  sections. 
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(b) 


(c) 


Figure  29  -  (a)  Far-field  radiation  patterns  vs.  wavelength  of  operation.  At  722  nm  (solid  blue  line),  714  nm  (dashed 
red  line)  and  690nm  (dotted  green  line),  directional  far-field  radiation  patterns  are  obtained,  pointing  at  20°,  18°,  and 
13°  respectively.  At  600nm,  the  guided- wave  mode  does  not  significantly  contribute  to  the  far-field  radiation,  (b) 
Calculated  three-dimensional  leaky-wave  radiation  pattern  at  the  wavelength  of  690nm.  (c)  Scanning  of  the  main 
lobe  radiation  pattern  (magnitude  and  main  direction)  versus  wavelength.  The  highlighted  region  corresponds  to 

leaky-wave  operation. 

Figure  29(a)  shows  the  corresponding  far-field  radiation  patterns  in  the  E  plane  at  various 
wavelengths.  It  is  seen  that,  in  the  leaky- wave  regime,  the  conical  beam  may  scan  the  angle  with 
frequency,  as  predicted  in  the  previous  sections.  The  patterns  show  a  significant  directivity  that 
may  be  tuned  by  changing  the  frequency  of  operation  (i.e.,  the  material  permittivity).  It  is  seen 
that,  consistent  with  the  previous  analytical  results,  better  directivity  is  achieved  for  radiation 
closer  to  the  chain  axis,  for  which  £  is  larger.  The  scan  of  the  main  lobe  direction  with 
frequency  confirms  the  forward  nature  of  these  longitudinal  leaky-wave  modes,  as  predicted  by 


55 


Approved  for  public  release;  distribution  unlimited. 


the  previous  analysis.  For  comparison,  the  radiation  at  A0  =600nm  is  very  poor,  due  to  the 

guided-wave  properties  of  the  chain  at  this  wavelength.  The  subdiffractive  nature  and 
subwavelength  period  of  the  chain  ensure  absence  of  significant  side  lobes.  These  results  confirm 
the  realistic  possibility  of  using  a  silver  nanoparticle  chain  as  a  leaky-wave  nanoantenna. 
Different  nanoparticle  size  and  geometry  may  be  used  to  tune  and  shift  the  leaky-wave  operation 
at  different  wavelengths. Figure  29b  reports  the  three-dimensional  far-field  radiation  pattern  at 
714  mn,  together  with  the  geometry  of  the  chain,  to  highlight  the  directive  conical  radiation  at  18 
degrees  from  the  chain  axis,  consistent  with  Fig.  29a.  Smaller  side  lobes  are  visible,  associated  to 
the  finite  length  of  the  chain.  As  reported  in  Fig.  29c,  the  nanoparticle  chain  supports  a  smooth 
linear  scanning  region  between  the  wavelengths  of  680  to  740  mn  (highlighted  in  the  figure), 
which  delimit  the  leaky-wave  operation  of  this  nanoantenna.  Tunability  and  beam  scanning  at  the 
same  frequency  may  be  envisioned  by  considering  electro-optical  materials  or  proper 
nonlinearities  in  the  nanoparticles. 

g.  Conclusions 

In  this  section,  we  have  provided  a  detailed  analysis  of  the  general  leaky-wave  radiation 
properties  of  linear  arrays  of  subwavelength  plasmonic  nanoparticles.  Using  closed-form 
analytical  dispersion  relations  for  real  and  complex  modes  supported  by  such  chain,  we  have 
analyzed  its  leaky-mode  properties  and  the  most  general  conditions  required  to  support  this 
regime  with  large  directivity  and  robust  frequency  response,  in  the  limit  of  subwavelength 
nanoparticles  composing  the  array. 

In  particular,  we  have  shown  that  the  longitudinal  polarization  is  the  best  candidate  for  achieving 
significantly  directive  conical  radiation  with  scanning  capabilities  in  this  regime,  and  that  the 
transverse  polarization  may  support  backward-wave  radiation.  We  have  also  considered  the 
effects  of  varying  the  center-to-center  distance,  the  material  properties  and  the  possible  presence 
of  material  dispersion  and  loss,  specializing  our  general  analysis  to  realistic  plasmonic  materials 
and  providing  comparison  with  analogous  dielectric  arrays.  Our  analysis  shows  that  plasmonic 
materials  may  provide  a  robust  route  to  leaky-wave  radiation  at  optical  frequencies,  adding  more 
flexibility  over  the  leaky-wave  properties  of  thin  plasmonic  films.  We  have  also  validated  our 
results  with  full-wave  simulations,  analyzing  the  leaky-wave  propagation  and  radiation  along 
silver  nanosphere  arrays  of  finite  length.  Our  full-wave  simulations  have  confirmed  that 
plasmonic  leaky-wave  nanoantennas  with  sub-diffractive  lateral  cross  section  may  indeed  lie 
within  the  realm  of  current  nanotechnology  and  may  be  applied  to  novel  devices  for  optical 
communications  and  computing. 


6.  Coupling  and  Guided  Propagation  along  Parallel  Arrays  of  Nanoparticles 

a.  Summary 

In  this  section,  we  derive  a  dynamic  closed-form  dispersion  relation  for  the  analysis  of  the  entire 
spectrum  of  guided  wave  propagation  along  coupled  parallel  linear  arrays  of  plasmonic 
nanoparticles,  operating  as  optical  “two-line”  waveguides.  Compared  to  linear  arrays  of 
nanoparticles,  our  results  suggest  that  these  waveguides  may  support  more  confined  beams  with 
comparable  or  even  longer  propagation  lengths,  operating  analogously  to  transmission-line 
segments  at  lower  frequencies.  Our  formulation  fully  takes  into  account  the  entire  dynamic 
interaction  among  the  infinite  number  of  nanoparticles  composing  the  parallel  arrays, 
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considering  also  realistic  presence  of  losses  and  the  frequency  dispersion  of  the  involved 
plasmonic  materials,  providing  physical  insights  into  the  guidance  properties  that  characterize 
this  geometry. 

b.  Introduction 

Linear  chains  of  plasmonic  (silver  or  gold)  nanoparticles  have  been  suggested  as  optical 
waveguides  in  several  recent  papers  [131]-[141],  as  outlined  in  the  previous  sections.  Owing  to 
design  flexibility  and  relatively  easy  construction  within  current  nanotechnology,  the  realization 
of  such  ultracompact  waveguides  has  been  thoroughly  studied  and  analyzed  in  the  past  few 
years.  However,  the  recent  experimental  realizations  of  such  devices  at  the  nanoscale  have 
revealed  challenges  due  to  severe  sensitivity  to  material  absorption  and  to  inherent  disorder.  The 
guided  beam  cannot  usually  travel  longer  than  few  nanoparticles  before  its  amplitude  is  lost  in 
the  noise.  This  is  mainly  due  to  the  fact  that  linear  arrays  of  small  nanoparticles  have  the 
property  to  concentrate  the  optical  energy  in  a  sub-wavelength  region  of  space,  in  large  part 
filled  by  lossy  metal.  If  this  is  indeed  appealing  in  tenns  of  power  concentration  and  for 
enhancing  the  nanoscale  interaction  with  light,  it  also  has  the  clear  disadvantage  of  strong 
sensitivity  to  material  and  radiation  losses.  In  general,  there  is  a  well-known  trade-off  between 
energy  concentration  and  confinement  and  robustness  to  loss  in  several  plasmonic  waveguide 
geometries  [142].  As  we  have  noted  in  [143],  a  bare  conducting  wire  at  low  frequencies  has 
analogous  limitations:  although  metals  are  much  more  conductive  and  less  lossy  in  radio 
frequencies,  connecting  two  points  in  a  regular  circuit  with  a  single  wire  would  still  produce 
unwanted  spurious  radiation  and  sensitivity  to  metal  absorption.  This  problem,  which  is  much 
amplified  at  optical  frequencies  due  to  the  poorer  conductivity  and  higher  loss  of  metals  in  the 
visible,  is  simply  approached  at  low  frequencies  by  closely  pairing  two  parallel  wires  (or,  which 
is  the  same,  placing  a  ground  plane  underneath  the  conducting  trace),  forming  the  well-known 
concept  of  a  transmission-line  that  provides  a  return  path  for  the  conduction  current. 
Analogously,  applying  the  nanocircuit  concepts  [144]-[145],  we  have  recently  put  forward  ideas 
to  realize  optical  nanotransmission-line  waveguides  in  different  geometries  [146] -[147],  which 
have  been  proven  to  be  more  robust  to  material  and  radiation  losses  and  may  provide  wider 
bandwidth  of  operation.  In  particular,  one  such  idea  consists  in  pairing  together  two  parallel 
arrays  of  plasmonic  nanoparticles,  suggesting  that  the  coupling  among  the  guided  modes  may 
improve  the  guidance  performance.  In  [143]  we  have  shown  that  this  is  indeed  the  case: 
operating  with  the  antisymmetric  longitudinal  mode,  such  parallel  chains  may  indeed  confine  the 
beam  in  the  background  region  between  the  chains,  leading  to  large  power  confinement  without 
significantly  affecting  the  robustness  to  material  absorption  and  radiation  losses  as  compared  to 
the  isolated  array  scenario.  In  particular,  we  have  shown  that  operating  with  these  modes  near  the 
light-line  would,  in  many  senses,  lead  to  operation  close  to  a  regular  transmission-line  at  low 
frequencies,  but  available  in  the  visible  regime. 

Here,  we  present  a  complete,  closed-form  dynamic  solution  for  the  dispersion  of  the  eigenmodes 
supported  by  parallel  chains,  fully  taking  into  account  the  coupling  among  the  infinite  number  of 
particles  composing  the  two-chain  array,  even  in  the  presence  of  material  absorption,  radiation 
losses  and  frequency  dispersion.  This  derivation  allows  us  to  discuss  the  complete  spectrum  of 
guided  modes  supported  by  this  geometry  and  analyze  the  differences  among  different 
polarizations  and  the  analogies  with  the  isolated  chain  geometry.  The  results  confirm  the  validity 
of  the  analogy  between  these  parallel  arrays  of  nanoparticles  and  optical  transmission  lines,  and 
they  provide  further  insights  into  the  operation  and  the  large  spectrum  of  modes  guided  by  these 
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paired  arrays  of  nanoparticles.  Applications  for  low-loss  optical  communications,  optical 
switching,  nonlinearity  enhancement  and  sub-wavelength  imaging  devices  are  envisioned. 


c.  Dispersion  Relations  for  Guided  Propagation 

Consider  the  geometry  of  Fig.  30,  i.e.,  two  identical  linear  arrays  of  plasmonic  nanoparticles  with 
radius  a ,  period  d  >  2a  and  interchain  distance  /  >  d  .  This  geometry  has  been  preliminarily 
analyzed  in  [143]  for  its  longitudinally  polarized  antisymmetric  guided  modes,  where  it  was 
shown  that  the  coupling  between  the  chains,  limited  in  that  analysis  to  its  dominant  contribution 
coming  from  the  averaged  current  density  on  the  chain  axes,  was  expected  to  generate  the 
splitting  of  the  regular  longitudinal  mode  guided  by  an  individual  linear  array  into  two  coexisting 
longitudinal  modes,  respectively,  with  symmetric  and  antisymmetric  field  distributions  and 
polarization  properties,  as  sketched  in  Fig.  30.  The  antisymmetric  mode  is  the  one  corresponding 
to  transmission-line  operation  [143],  as  outlined  in  the  introduction,  for  which  two  antiparallel 
displacement  current  flows  are  supported  by  the  parallel  chains.  A  similar  modal  propagation  has 
been  analyzed  in  [137]  for  a  related  distinct  geometry,  consisting  of  longitudinal  dipoles  placed 
over  a  perfectly  conducting  plane.  Also  our  analysis  of  quadrupolar  chains  [147]  may,  in  the 
limit  of  /  — >  0 ,  have  some  analogies  with  this  antisymmetric  operation.  In  the  following,  we 
rigorously  approach  the  general  problem  of  modal  dispersion  along  the  parallel  chains  of  Fig.  30, 
extending  our  general  analysis  in  [140]  valid  for  one  isolated  chain.  Our  formulation  fully  takes 
into  account  the  entire  coupling  among  the  infinite  nanoparticles  composing  the  pair  of  arrays 
and  the  possible  presence  of  material  absorption,  radiation  losses  and  frequency  dispersion. 
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Figure  30  -  Geometry  of  the  problem:  a  pair  of  linear  arrays  of  plasmonic  nanoparticles  as  an  optical  two-line 
waveguide.  The  sketch  reports  the  polarization  properties  of  quasi-longitudinal  modes  with  symmetric  and 

antisymmetric  properties. 

We  model  each  nanoparticle  in  the  coupled  array  of  Fig.  30  as  a  polarizable  dipole  with 
polarizability  a  ,  an  assumption  that  is  valid  as  long  as  a  «  Xb ,  with  Xh  being  the  wavelength  of 

operation  in  the  background  material,  consistent  with  the  assumptions  in  the  previous  sections. 
For  simplicity,  we  assume  a  scalar  polarizability,  implying  that  the  particles  are  isotropic 
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(nanospheres,  easy  to  realize  as  colloidal  metal  particles),  or  for  more  general  shapes  focusing  on 
one  specific  field  polarization.  It  is  relevant  to  stress  that  the  dipolar  approximation  represents  a 
good  assumption  for  small  nanoparticles,  and  in  particular  in  the  case  of  nanospheres,  due  to 
their  inherent  symmetries.  We  have  verified  that  this  approximation  holds  very  well  even  in  the 
limit  of  very  small  gaps,  as  those  considered  in  the  following  examples  [140].  The  use  of 
additional  multipolar  orders,  as  considered,  e.g.,  in  [148],  may  increase  the  accuracy  of  the 
calculation,  but  also  complicate  some  of  the  physical  insights  outlined  in  the  following. 

For  a  single  isolated  chain,  the  spectrum  of  supported  eigenmodes  may  be  split  into  longitudinal 
and  transverse  polarization  with  respect  to  the  chain  axis  x  [140]  (see  previous  sections).  In 
particular,  for  e‘/lx  propagation,  the  corresponding  guided  wave  number  [3  satisfies  the 
following  closed-fonn  dispersion  relations,  respectively,  for  longitudinal  and  transverse  modes: 


L=3d 


-a  1  =  0 


T=  —  d 
2 


2  /  —  '  f  /■  /  7  \  •  i  r  (  n  i\  l 


fi(P’d)~id  f2{P’d)-d2fl{/3,d) 


-a  1  =  0 


(34) 


where  fN[(3,  d^  =  LiN  j  +  LiN  (e  ^  j ,  LiN  (z)  =  is  the  polylogarithm  function 

\  \  /  k=]  k 

of  order  N  [149]  and  all  the  quantities  have  been  normalized,  consistent  with  [140],  as  d  =  kbd  , 
(3  =  P  /  kb,  a  -  k\a  /  ( 6  neb ) ,  with  kh  =  2 jz  I  Ah  being  the  background  wave  number  and  sb  the 
corresponding  permittivity. 

These  equations  fully  take  into  account  the  dynamic  coupling  among  the  infinite  number  of 
particles  composing  the  linear  chain.  They  are  real-valued  for  lossless  particles  (for  which 
Im  [a  1  ]  =  — 1  [140]),  supporting  guided  modes  with  (3  >  1,  but  they  are  also  fully  valid  in  the 

complex  domain  when  realistic  material  losses  are  considered,  making  it  possible  to  evaluate  the 
realistic  damping  factors  associated  with  material  absorption  and  radiation  losses.  They  can  be 
applied  also  to  the  leaky- wave  modal  regime,  for  which  Re  [/?  J  <  1  and  the  chain  radiates  as  an 
antenna  in  the  background  region  [150]. 

When  /  is  finite  in  Fig.  30,  i.e.,  there  are  two  parallel  chains,  their  mutual  coupling  produces  a 
modification  of  their  guidance  properties,  which  may  be  taken  into  account  by  considering  the 
polarization  fields  induced  by  the  electric  field  from  each  chain  on  the  other.  The  fields  radiated 
by  each  chain  may  be  expanded  into  cylindrical  waves,  allowing  us  to  write  the  general  closed- 
form  expressions  for  the  coupling  coefficients  between  the  two  chains. 

Without  loss  of  generality,  we  can  assume  that  the  particles  composing  the  first  chain,  located  at 
y  =  0,  are  polarized  by  an  eigenmodal  wave  with  dipole  moments  p,e'/w  ,  where  m  e  3  is  the 
integer  index  for  each  nanoparticle  of  the  chain.  The  equivalent  current  distribution  on  the  x  axis 
may  be  written  as: 

oo 

J  (x)  =  -ifflpj  X  elpmd5(x-md)  ,  (35) 
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where  S  (.)  is  the  Dirac  delta  function.  The  fields  radiated  by  such  current  distribution  may  be 

expanded  into  cylindrical  waves  and  may  be  used  to  evaluate  the  coupling  coefficients  between 
one  chain  and  the  other,  with  dipole  moments  p,e,/w  located  at  y  =  l ,  yielding: 


O  oo 

C  =-=  V  b2K,Jb  H 

xx  i  m  0  [_  m  J 
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where  bm  = 

l  =kbl .  The  generic  coupling  coefficient  C  expresses  the  polarization  along  j  on  one  chain 

induced  by  the  /-polarized  dipoles  on  the  other  chain.  The  summations  in  (36)  have  very  fast 
convergence,  and  the  dominant  term  ( m  —  0 )  is  usually  sufficient  to  take  into  account  the 
dominant  contribution  to  the  coupling,  an  approximation  that  is  consistent  with  the  approach  we 
used  in  [143].  The  numerical  results  shown  in  the  following  sections  have  been  obtained  by 
considering  the  first  ten  terms  in  the  summations  (36),  even  though  full  convergence  has  been 
usually  achieved  after  the  first  one  or  two  terms.  The  form  of  the  coupling  coefficients  in  (36) 
ensures  that  longitudinal  (directed  along  x)  and  transverse  modes  polarized  along  y  are  coupled 
through  C  ,  whereas  transverse  modes  polarized  along  z  are  not  coupled  with  the  orthogonal 
polarizations. 
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-1  and  Km[]  are  the  modified  Bessel  functions  of  order  m  and 


The  complete  closed-fonn  dispersion  relation  for  the  eigenmodes  supported  by  the  parallel 
chains  may  be  written  as: 
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or,  in  a  more  compact  form: 

[(z  ±  cv* )  (r + c„, ) + cA>  2  ]  (r  ±  czz )  =  o . 


(37) 


(38) 
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Frequency  [THz] 


Figure  31  -  Modal  dispersion  for  the  quasi-longitudinal  modes  supported  by  two  parallel  chains  with  interchain 
distance  l  =  50  nm  .  The  dispersions  are  compared  to  that  of  an  isolated  chain  (thin  solid  line).  The  lighter  blue 
shadowed  region  refers  to  leaky-wave  propagation,  whereas  the  darker  orange  region  refers  to  the  first  Bragg  stop- 

band. 
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Figure  32  -  Similar  as  in  Fig.  31,  modal  dispersion  for  the  quasi-longitudinal  modes  supported  by  two  parallel 

chains  with  interchain  distance  l  =  30  nm  . 

The  left-hand  side  in  Eq.  (38)  consists  of  the  product  of  two  terms:  the  first  determines  the 
dispersion  of  the  coupled  modes  polarized  in  the  xy  plane  (among  which  the  quasi-longitudinal 
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antisymmetric  modes  considered  in  [143]),  whereas  the  second  determines  the  purely-transverse 
modes  polarized  along  z  .  This  dispersion  equation  is  completely  general  and  it  fully  takes  into 
account  the  entire  dynamic  interaction  among  the  infinite  particles  composing  the  two  parallel 
chains.  Since  the  coupling  coefficients  (36)  tend  rapidly  to  zero  for  increased  /,  Eq.  (38) 
represents  the  perturbation  of  the  original  transverse  and  longitudinal  modes  supported  by  the 
two  linear  chains  independently  given  by  L  =  0  and  T  —  0  respectively  [140],  produced  by  the 
coupling  coefficients  C .  In  particular,  it  is  seen  that  each  of  the  three  orthogonal  polarizations 
(along  x,  y,  z)  splits  into  two  branches  due  to  the  coupling  between  the  chains,  one  with 
symmetric  and  the  other  with  antisymmetric  properties  (consistent  with  the  sketch  in  Fig.  30), 
leading  to  six  modal  branches  of  guided  modes,  some  of  which  supported  at  the  same  frequency. 
In  particular,  the  modes  in  the  xy  plane  are  mixed  together  (i.e.,  the  parallel  chains  do  not 
support  purely  longitudinal  or  purely  y  -  polarized  modes,  but  they  do  support  purely  transverse 
z- polarized  modes).  In  the  limit  of  lossless  particles,  since  L  and  T  are  real  for  any  ft  >  1 
[140],  by  inspecting  Eq.  (38)  we  notice  that  the  parallel  chains  still  support  lossless  guided 
propagation  for  any  \<  fd  <n  !  d  An  the  following,  we  analyze  in  details  the  modal  properties  of 
this  setup  in  its  different  regimes  of  operation. 


d.  Guided  Modes  of  Parallel  Chains  of  Silver  Nanospheres 

In  this  section  we  consider  the  different  regimes  of  guided  propagation  supported  by  the  parallel 
chains  of  Fig.  30,  considering  realistic  optical  materials  composing  the  plasmonic  nanoparticles. 
In  the  case  of  a  chain  of  homogeneous  spherical  particles  of  radius  a  and  permittivity 
e  =  er+  i  si ,  their  normalized  polarizability  satisfies  the  following  relations  [140]: 


Re[a  '] 


er-£b 

£b{hay 
\2  2 
£r~£b)  +  £i 


(39) 


This  polarizability  model  is  analogous  to  the  classic  quasi-static  polarizability  definition  for  a 
small  sphere,  with  the  addition  of  radiation  loss  (as  from  the  -1  term  in  Im^a_1J),  to  comply 

with  energy  conservation.  Since  the  guided  modes  are  perturbations  of  the  longitudinal  and 
transverse  modes  supported  by  the  isolated  chains,  there  is  no  need  to  analyze  here  again  in  full 
detail  how  variations  in  the  chain  geometry,  i.e.,  in  a,  d  and/or  the  involved  materials,  may 
affect  the  guidance  of  the  parallel  chains,  since  we  have  already  extensively  studied  how  these 
changes  affect  the  guidance  of  isolated  chains  in  [140].  In  the  following,  therefore,  we  focus  on 
one  specific  realistic  design  of  the  chains  and  we  employ  the  exact  fonnulation  developed  in  the 
previous  section  to  characterize  the  modal  properties  of  two  of  such  parallel  arrays  coupled 
together.  In  the  following,  we  focus  on  colloidal  silver  nanospheres  embedded  in  a  glass 
background  ( sh  =  2.38  sQ  ).We  use  experimental  data  available  in  the  literature  to  model  the  silver 
pennittivity  at  optical  frequencies  [27]  and  we  assume  a  =  IQnm  and  <7  =  2 1  n m  for  the  two 
chains.  In  general,  we  can  predict  that  larger  separation  distances  for  same  particle  size  are 
expected  to  weaken  the  guidance  properties  of  the  array  and  reduce  the  bandwidth  of  operation 
and  robustness  to  loss,  whereas  particles  with  larger  size  have  the  opposite  effect. 


63 


Approved  for  public  release;  distribution  unlimited. 


i.  Quasi-longitudinal  propagation  (forward  modes) 

An  isolated  linear  chain  of  plasmonic  nanoparticles  supports  forward-wave  longitudinal  guided 
modes  (x- polarized),  satisfying  the  dispersion  relation  L-  0,  over  the  frequency  regime  for 
which: 

6[c/3(J+;t)  +  dCl2[d  +  7r)\<  <73  Re[a_1]  <  3[|(3)  +  C/3  (2j)  +  d  Cl2[ld)\,  (40) 

where  ClN(d)  are  Clausen’s  functions  [149]  and  £,(.)  is  the  Riemann  Zeta  function.  For  the 
case  at  hand  (silver  nanoparticles,  a  =  10nm  and  d  =  21nm),  such  modal  regime  is  supported 
over  a  relatively  wide  range  of  frequencies  between  550  THz  and  850  THz  ,  as  shown  in  Fig.  31 
(thin  solid  black  line).  In  particular,  in  the  figure  we  plot:  (a)  the  real  and  (b)  the  imaginary  parts 
of  the  normalized  (5  and  (c)  the  propagation  length,  i.e.,  the  distance  traveled  by  the  guided 

mode  before  its  amplitude  is  e~l  of  the  original  value,  which  is  equal  to  Im  [/?]  ' .  The  shadowed 
regions  at  the  sides  of  the  plots  delimit  the  leaky-wave  region  (left-side,  lighter  blue  shadowed 
region),  for  which  Re[/?]<1  and  the  mode  radiates  in  the  background  region  [154],  and  the 

stop-band  region  (right  side,  darker  shadow,  brown),  where  Re[/?]  =  ;r/(i  when  lossless 
particles  are  considered  and  the  mode  is  evanescent  in  nature.  In  between  these  two  regions,  as 
defined  by  Eq.  (40),  the  modes  are  guided  and  Im[/?],  i.e.,  the  damping  factor,  is  only 
associated  with  material  losses,  since  in  the  limit  of  lossless  particles  the  mode  would  not  radiate 
and  kn[/?]  =  0.  In  the  leaky- wave  region  (lighter  blue  shadow)  the  damping  is  larger,  due  to 

additional  radiation  losses  [150],  whereas  in  the  stop-band  (darker  shadow)  the  mode  does  not 
propagate  and  it  is  reflected  back  by  the  chain  due  to  Bragg  reflection  (the  ideal  Bragg  stop-band 
line  R e[/?]  =  ;r/c/  is  shown  in  Fig.  31  as  the  thin  pink  line).  Near  the  light  line  (Re[/?]  =  1)  the 

mode  is  poorly  guided  by  an  isolated  chain,  but  its  propagation  length  may  reach  relatively  large 
values,  around  1  pm  . 

In  the  same  figure  we  plot  the  variation  of  these  dispersion  diagrams  in  the  case  of  two  coupled 
parallel  chains  separated  by  a  finite  distance  / .  In  this  case,  the  longitudinal  modes  are  coupled 
with  each  other,  also  polarizing  the  chains  with  a  small  transverse  component  along  y , 
consistent  with  the  value  of  C  .  The  longitudinal  mode  dispersion  splits  into  two  quasi¬ 
longitudinal  branches,  one  with  symmetric  and  the  other  with  antisymmetric  properties  with 
respect  to  x-  polarization  (as  sketched  in  Fig.  30).  The  two  modes  satisfy,  respectively,  the 
following  dispersion  relations,  consistent  with  Eq.(38): 

sym  •  (i+cj(r-cj+cy=0 

antisym :  (I-C„)(r  +  C„  ,)+Cj  =0 

providing  the  following  constraints  on  the  polarization  eigenvectors  for  the  two  chains  [obtained 
by  calculating  the  eigen-vectors  associated  with  Eq.  (37)]: 
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(42) 


sym : 


antisym : 


|prx  =  p2-x 

|pry  =  -p2-y 

rPl-x=-p2-x 

lpry  =  p2-y 


Figure  33  -  (Color  online).  Magnetic  field  distribution  (snapshot  in  time)  for  the  chains  of  Fig.  32  at  frequency 
/  =  585  TIIz  .  (a)  antisymmetric  mode,  (b)  symmetric  mode,  (c)  isolated  chain.  All  the  plots  are  drawn  with  the 
same  color  scale  bar  (normalized  to  the  modal  amplitude  at  the  left  of  the  figure).  The  total  length  of  the  simulated 

region  is  21b . 


Figure  31  shows  as  a  first  example  the  dispersion  of  symmetric  and  antisymmetric  modes  for 
/  =  50  nm  .  It  is  noticed  that  the  small  coupling  between  the  chains  slightly  perturbs  the  dispersion 
of  the  modes,  causing  the  antisymmetric  mode  (blue  dashed  line,  with  polarization  currents 
oppositely  flowing  along  the  chains)  to  have  slightly  larger  real  and  imaginary  parts  of  /?  with 
respect  to  the  unperturbed  longitudinal  mode  supported  by  an  isolated  chain  (light  solid  line)  in 
the  guided  regime.  Conversely,  the  symmetric  mode  (thick  red  solid  line)  supports  slightly  lower 
values  of  Re  J  .  The  perturbation  is  stronger  near  the  light  line  and  in  the  leaky- wave  region, 

since  the  mode  is  less  confined  around  each  chain  in  this  regime.  The  symmetric  operation 
allows  an  increase  of  the  propagation  length  of  up  to  1.5  /um ,  since  the  coupling  between  the 
parallel  chains  with  polarization  currents  flowing  in  the  same  direction  can  boost  up  the  mode. 
On  the  other  hand,  the  antisymmetric  operation  has  slightly  lower  propagation  lengths,  but  this  is 
accompanied  with  the  important  advantage  of  much  stronger  field  confinement,  as  we  note  in  the 
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following.  The  derivative  a  Re[/?]  /  dco  >  0  ensures  that  the  modes  supported  in  this  regime  are 
all  forward- wave,  and  this  is  also  confirmed  by  the  condition  Im[/?]  >  0,  which  ensures  that 
phase  and  group  velocity  are  parallel  with  each  other  for  both  modes. 

As  an  aside,  it  should  be  noted  that  in  the  leaky-wave  region  (blue  lighter  shadow  in  the  left)  the 
forward- wave  modes  are  improper  in  nature  [154],  implying  that  the  dominant  cylindrical  wave 
radiated  by  the  chain  grows  with  the  distance  from  the  chain  instead  of  decaying.  This  implies 
that  for  a  correct  evaluation  of  the  modal  properties  and  the  field  distribution  generated  in  this 
forward-leaky  mode  regime,  the  formulas  of  Eq.  (36)  for  the  index  m  —  0  need  to  be  corrected 
using  the  Hankel  functions  of  second  order  instead  of  the  modified  Kn  functions. 

Figure  32  shows  analogous  results  for  closer  chains,  with  l  =  30  nm.  It  is  seen  that  the 
perturbation  from  the  isolated  chain  is  now  stronger  and  the  coupling  between  the  modes 
generates  some  isolated  resonant  regions  of  stronger  absorption,  which  are  associated  with 
stronger  transversely  polarized  components  of  the  field.  Still,  near  the  light  line  propagation 
lengths  are  relatively  large. 

Figure  33  shows  the  calculated  orthogonal  magnetic  field  distribution  (snapshot  in  time)  on  the 
xy  plane  for  the  modes  supported  by  the  chains  of  Fig.  32  (/  =  30 nm )  at  the  frequency 
/  =  585  THz ,  near  the  light  line.  The  figure  emphasizes  how  the  modal  distribution  is  quite 
different  in  the  three  scenarios,  even  if  the  guided  wave  numbers  are  similar.  Fig.  33a 
corresponds  to  antisymmetric  propagation,  for  which  the  two  chains  support  the  eigenvector 
polarizations  p,  =  x  +  (0.14z -0.008)y  ,  p2  =  -x  + (0.14/ -0.008) y ,  consistent  with  Eq.  (42). 

The  corresponding  nonnalized  wave  number  at  this  frequency  is  fi  =1.38  +  / 0.1 .  It  can  be 

seen  how  the  magnetic  field  is  very  much  confined  in  the  tiny  background  region  delimited  by 
the  two  chains,  similar  to  the  field  propagation  in  a  regular  transmission-line  at  low  frequencies. 
Also  the  electric  field  is  mainly  transverse  in  the  region  between  the  chains,  supporting  the 
transverse  electromagnetic  configuration,  again  typical  of  a  transmission-line  mode.  This  regime 
of  operation,  whose  interesting  properties  we  have  already  described  in  details  in  [143],  may  lead 
to  relatively  low-loss  optical  guidance  with  ultra-confined  properties  in  the  space  between  the 
chains,  similar  to  an  optical  nano  transmission- line.  This  functionality  may  be  particularly 
appealing  for  optical  switching,  nanoscale  field  interaction,  sensing  and  nonlinearity 
enhancements. 

Fig.  33b,  on  the  other  hand,  refers  to  the  symmetric  mode  for  the  same  parallel  chains.  In  this 
case  P[  =  x  +  (0.08/-0.003)y  ,  p2  =  x-(0.08/-0.003)y  and  j3sym  =1.13  +  / 0.053  .  The  currents 

flowing  along  the  chains  are  now  parallel  with  each  other,  producing  fields  very  much  spread  all 
around  the  outside  background  region  and  weak  field  concentration  in  between  them.  This 
operation  is  equivalent  to  two  parallel  current  flows,  leading  to  small  fields  in  between  them.  A 
single  linear  chain  has  analogous  guidance  properties,  shown  in  Fig.  33c  (for  comparison,  in  this 
third  example  the  chain  is  positioned  at  the  same  location  as  the  lower  chain  in  the  other  two 
panels).  In  this  case  the  mode  is  purely  longitudinal  and  /?singlc  =1.18  +  / 0.072 ,  implying  weak 
guidance. 

Comparing  the  three  field  plots  (notice  that  for  fair  comparison  they  have  been  calculated  with 
the  same  color  scale  and  under  the  same  initial  amplitude  excitation),  it  becomes  evident  that  the 
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antisymmetric  longitudinal  operation  allows  for  a  much  stronger  confinement  of  the  field,  with 
comparable  propagation  length.  By  field  confinement  here,  we  consider  the  field  decay  away 
from  the  array  axis,  which  is  always  of  a  higher  rate  for  the  anti-symmetric  longitudinal  mode 
compared  to  an  isolated  chain  or  the  symmetric  mode.  By  increasing  the  distance  between  the 
chains,  as  in  the  examples  of  Fig.  31,  we  achieve  similar  confinement  in  the  region  between  the 
chains  with  reduced  attenuation.  In  addition  to  a  quantitatively  larger  field  confinement  achieved 
with  the  antisymmetric  mode  (compared  to  the  isolated  chain,  despite  the  larger  transverse 
physical  cross  section  of  the  parallel  chain  geometry),  this  operation  ensures  larger  confinement 
in  between  the  two  chains,  which  is  particularly  appealing  for  applications  aiming  at  enhancing 
the  nanoscale  optical  interaction. 

These  properties  are  not  only  limited  to  the  modes  operating  near  the  light  line,  but  they  are  also 
valid  for  higher  frequencies  and  more  confined  modes.  For  instance,  in  Fig.  34  we  show  the 
magnetic  field  plots  for  the  same  chains,  operating  at  /  =  680  THz .  At  these  frequencies,  as  seen 
in  Fig.  32,  the  three  cases  have  similar  levels  of  absorption  and  more  confined  slow-wave  modes. 
The  antisymmetric  excitation  is  characterized  in  this  case  by  =x  +  (0.37  /'-  0.038) y, 

p2  =  -X  +  (0.37/  -0.038) y  and  J3asym  =  2.42  +  i 0.091 .  Its  field  distribution  (Fig.  34a)  still  shows 

strong  confinement  between  the  two  chains,  where  a  “quasi-uniform”  magnetic  field  may 
propagate  as  if  guided  by  a  transmission-line.  The  wave  is  slower  than  in  the  case  of  Fig.  4,  due 
to  increased  Re  [/?]  ,  but  the  level  of  absorption  is  still  quite  good  and  the  mode  can  propagate 

for  over  two  wavelengths  with  no  strong  attenuation.  The  symmetric  operation,  for  which 
Pj  =  x-(0.158/' -0.082)y  ,  p2  =  x  +  (0.158/'-0.082)y  and  J3sym  =1.9  +  i 0.074,  once  again 

provides  worse  field  confinement,  as  expected.  In  this  case  (Fig.  34b)  the  field  is  spread  around 
the  chains  and  is  very  weak  in  the  region  between  the  two  chains.  Similar  spreading  is  noticeable 
in  the  single  isolated  chain  configuration  of  Fig.  34c,  with  /?single  =  2.06  +  / 0.077  .  We  note  that 

the  field  spreading  in  the  region  around  the  chains  would  also  be  more  sensitive  to  radiation 
losses  produced  by  disorder  and  technological  imperfections.  This  is  consistent  with  our 
findings  in  [141],  in  which  we  have  quantitatively  modeled  the  possible  presence  of  disorder 
along  such  periodic  arrays  in  tenns  of  effective  additional  loss  in  the  polarizability  coefficient. 
We  predict  that  the  antisymmetric  transmission-line  operation  of  the  parallel  chains  may  produce 
more  robust  optical  guidance  confined  in  the  region  between  the  chains  compared  to  the 
symmetric  operation  or  the  isolated  chain. 

From  the  previous  examples,  it  is  evident  that  in  this  regime  the  modes  guided  by  the  parallel 
chains  are  quasi-longitudinal  with  a  spurious  transverse  polarization,  arising  from  the  coupling, 
which  is  nearly  90°  out  of  phase  with  respect  to  the  longitudinal  polarization.  In  Fig.  35,  for  the 
parallel  chains  of  Fig.  31  and  32  we  have  calculated  the  level  of  transverse  cross-polarization 
(defined  as  the  ratio  of  the  longitudinal  to  transverse  component  of  p  )  induced  on  the  particles 
due  to  coupling,  as  a  function  of  frequency.  It  is  evident  that  its  level  increases  for  closer  chains, 
as  expected,  and  it  is  larger  for  antisymmetric  modes.  In  the  region  of  enhanced  absorption  that 
we  have  noticed  in  Fig.  32,  the  corresponding  level  of  cross-polarization  is  also  very  high, 
implying  that  at  some  resonance  frequencies  the  transverse  polarization  may  be  even  higher  than 
the  longitudinal  one,  noticeably  affecting  the  chain  guidance.  As  expected,  in  these  regimes  the 
losses  are  inherently  larger.  The  coupling  is  minimal  near  the  light  line  and  in  the  leaky-wave 
and  stop-band  regimes,  and  reaches  its  maximum  somewhere  inside  the  guidance  region,  whose 
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position  in  frequency  varies  depending  on  the  distance  between  the  chains  and  the  mode  of 
operation. 


Figure  34  -  Similar  to  Fig.  33,  magnetic  field  distribution  (snapshot  in  time)  for  the  chains  of  Fig.  32  at  frequency 
/  =  680  THz  :  (a)  antisymmetric  mode,  (b)  symmetric  mode,  (c)  isolated  chain. 


- Symmetric  I  =  30  nm 

- Antisymmetric  I  =  30  nm 

- -  Symmetric  I  =  50  nm 


Figure  35  -Magnitude  of  the  transverse  cross-polarization  for  the  chains  of  Figs.  31  and  32,  operating  in  their  quasi¬ 
longitudinal  forward-wave  regime. 


68 


Approved  for  public  release;  distribution  unlimited. 


ii.  Quasi-transverse  y  -polarized propagation  (backward  modes) 


We  have  discussed  in  [140]  that  a  single  isolated  linear  chain  may  also  support  transversely 
polarized  guided  modes,  satisfying  the  exact  dispersion  relation  T  =  0 .  In  this  case,  the  condition 
on  the  particle  polarizability  is: 


d  3a!  <  d 3  Re  \  a  1 1  <  -3 

min 


C/3(u?  +  7r^j  +  d  CI2(d  +  n'j- d2 C/j  (g?  +  7r^j  , 


(43) 


where  ajn  is  defined  in  [140],  In  this  regime  the  isolated  chain  always  supports  two  distinct 

modes  at  the  same  frequency,  both  with  the  same  transverse  polarization:  one  is  guided  along  the 
chain  and  has  backward-wave  properties,  the  other  is  weakly  guided,  with  forward-wave 
properties  and  Re  [/? ]  - 1  (this  eigenmode  is  basically  a  simple  plane  wave  traveling  in  the 

background  region,  weakly  polarizing  the  nanoparticles.  This  is  not  of  interest  for  guidance 
purposes  [140],  but  it  is  still  discussed  here  for  sake  of  completeness).  For  the  geometry  at  hand, 
transversely-polarized  propagation  is  supported  over  the  frequencies  between  650  THz  and 
800  THz  ,  in  part  overlapping  with  the  longitudinally-polarized  regime,  as  shown  in  Fig.  36  (thin 
black  line),  consistent  with  Eq.  (43).  Remarkable  differences  are  noticed  between  longitudinal 
and  transverse  polarization:  the  confined  transverse  mode  is  backward  in  nature,  explaining  the 
negative  slope  of  Re[/?]  versus  frequency  and  the  negative  sign  of  Im[/?].  Correspondingly, 

the  bandwidth  is  more  limited  and  losses  are  higher,  as  is  usually  the  case  for  backward-wave 
waveguides.  As  a  consequence,  the  leaky- wave  operation  arises  now  at  the  upper  boundary  of 
the  guided  regime  (lighter  shadowed  region  in  Fig.  36)  [150],  whereas  the  Bragg  stop-band  is 
positioned  at  the  lower-end  of  the  guidance  band  (darker  shadow). 

Due  to  modal  coupling  in  the  xy  plane,  when  the  coupling  between  parallel  chains  is  considered, 
the  quasi-transverse  modes  still  satisfy  the  dispersion  relations  (41)  and  the  polarization 
eigenvectors  obey  the  same  relations  (42).  It  should  be  noticed,  however,  that  in  this  regime  the 
modes  are  quasi-transverse,  and  therefore  the  antisymmetric  mode  now  corresponds  to  parallel 
y-  polarized  chains,  whereas  the  symmetric  mode  supports  anti-parallel  polarization  along  y  , 
consistent  with  (42),  as  sketched  in  Fig.  34. 

Figure  36  shows  the  dispersion  of  symmetric  and  antisymmetric  quasi-transverse  modes  for 
l  =  50  nm.  Once  again,  the  relatively  small  coupling  between  the  chains  produces  a  minor 
perturbation  of  the  original  backward-wave  purely-transverse  mode,  which  causes  the 
antisymmetric  mode  to  have  slightly  lower  real  and  slightly  larger  imaginary  part  of  J3 . 

Conversely,  the  symmetric  mode  supports  slightly  larger  values  of  Re[/?]  .  As  in  the  previous 

section,  the  coupling  is  stronger  near  the  light  line  and  in  the  leaky-wave  region  (lighter  shadow), 
as  expected.  Also  in  this  scenario  the  symmetric  operation  allows  longer  propagation  lengths, 
even  if  in  this  case  the  y  -  polarized  currents  are  oppositely  oriented  (see  Fig.  37). 

Compared  to  quasi-longitudinal  forward  modes,  the  propagation  length  is  significantly  reduced 
for  these  confined  modes,  due  to  the  backward-wave  nature  of  these  modes  and  their  transverse 
polarization.  Of  course,  following  the  results  in  [140],  the  propagation  length  may  be  somewhat 
increased  and  optimized  by  increasing  the  size  of  the  nanoparticles  and/or  reducing  the 
interparticle  distance  d  .  Compared  to  backward-mode  supported  by  isolated  chains,  these  results 
show  that  the  coupling  between  two  parallel  chains  may  increase  the  propagation  length  of 
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backward-wave  optical  nanowaveguides,  simultaneously  improving  their  field  confinement  in 
the  space  between  the  two  lines,  which  may  be  of  interest  for  several  applications. 
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Figure  36  -  Analogous  to  Fig.  31,  modal  dispersion  for  the  quasi-  transverse  y  -  polarized  modes  supported  by  two 
parallel  chains  with  interchain  distance  l  =  50  nm  .  In  the  antisymmetric  polarization,  as  well  as  in  the  isolated  chain, 
it  is  evident  the  presence  of  a  second  quasi-transverse  mode  with  weakly  guided  properties,  explaining  the  presence 

of  two  distinct  branches. 
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Figure  38  shows  analogous  results  in  the  case  of  closer  chains  (/  =  30  nm  ).  Also  in  this  case,  the 
perturbation  from  the  isolated  chain  is  stronger  and  the  bandwidth  of  backward  operation  may  be 
increased  by  using  two  parallel  chains  in  the  symmetric  mode.  Flere  leaky  modes  (lighter 
shadow)  are  proper  in  nature  and  therefore  Eqs.  (36)  also  apply  to  this  regime  in  the  way  they  are 
written.  Both  in  Fig.  36  and  38,  for  completeness,  we  also  show  the  modal  branch  associated 
with  the  weakly  guided  forward-wave  transverse  mode,  which  is  visible  very  close  to  the  light 
line  (thin  solid  line  very  close  to  the  Re[/?]  =  1  line  in  both  plots).  Consistent  with  its  forward- 

wave  properties,  Im[/?]  >  0  for  this  mode.  As  outlined  above,  this  mode  is  of  minor  interest  for 

guidance  purposes,  since  it  is  a  minor  perturbation  of  a  plane  wave  traveling  in  the  background 
region,  very  weakly  affected  by  the  presence  of  the  chains.  It  is  noticed,  as  expected,  that  this 
second  branch  is  present  only  for  the  antisymmetric  modes,  whose  y  -  polarization  is  in  the  same 
direction  for  both  chains.  Indeed,  in  Figs.  36  and  38  one  can  see  a  second  branch  for  the  blue 
dashed  lines,  near  the  light  line,  corresponding  to  the  weakly  guided  forward  mode  traveling  in 
the  background  (not  of  interest  here). 

Figure  39  shows  the  magnetic  field  for  these  backward-wave  modes  as  in  Fig.  38  at  the 
frequency  /  =  100THz .  In  the  antisymmetric  case  (Fig.  39a)  Pj  =  (0.076-0. 12/)x  +  y , 

p2  = -(0.076-0. 12/)x  +  y  and  J3  =3.95-/0.46;  in  the  symmetric  case 
p,  =(0.02-0. 016/)x-y,  p2  =(0.02-0. 016z)x  +  y  and  fisym  =  5.88-/0.72 ;  for  the  isolated 
chain  /?sing]e  =  5.036-/0.42 .  The  field  distributions  in  some  senses  resemble  the  one  for  quasi¬ 
longitudinal  modes,  but  the  presence  of  a  dominant  transverse  polarization  does  not  allow  an 
analogous  strong  transmission-line  confinement  in  this  backward-wave  regime  for  the 
antisymmetric  modes.  Still,  the  plots  confirm  that  relatively  long  backward-wave  propagation 
(over  one  wavelength)  is  achievable  using  coupled  parallel  chains. 

Figure  40  shows  the  level  of  longitudinal  cross-polarization  for  the  chains  of  Figs.  36  and  38.  In 
this  scenario,  the  cross-polarization  is  in  general  lower  than  for  quasi-longitudinal  modes  and  it 
is  stronger  for  symmetric  modes.  Once  again,  the  cross-polarization  is  stronger  for  closely 
coupled  chains  and  it  has  some  resonant  peaks  in  the  middle  of  the  guidance  region,  for  which 
the  damping  is  correspondingly  increased. 
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Figure  37  -  Polarization  properties  of  the  symmetric  and  antisymmetric  quasi-transverse  modes  supported  by 

parallel  chains  as  in  Fig.  30. 
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Figure  38  -  Similar  to  Fig.  36,  modal  dispersion  for  the  quasi-  transverse  y  -  polarized  modes  supported  by  two 

parallel  chains  with  interchain  distance  /  =  30  nm  . 

Hi.  Purely  transverse  z-polarized  propagation  (backward  modes) 

When  the  chains  are  polarized  along  z  the  supported  modes  are  purely  transverse,  consistent 
with  (37).  Due  to  symmetry,  the  properties  for  isolated  chains  are  identical  to  those  described  in 
the  previous  section,  and  therefore  here  we  discuss  how  the  coupling  may  have  a  different  effect 
on  the  backward-wave  guidance  properties  in  this  polarization.  The  coupling  coefficient  Czz 
splits  the  transverse  modal  branch  of  propagation  into  two  modes,  with  dispersion  relations: 
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sym\  T  +  CZ_  =  0 

antisym  :  T  -  C2Z  =  0 ' 


(44) 


providing  the  following  constraints  on  the  polarization  eigenvectors  for  the  two  chains: 
sym:  p,z  =  p2z 

antisym  :  p  z  =  -p2  •  z 


(45) 


Figure  39  -  Magnetic  field  distribution  (snapshot  in  time)  for  the  chains  of  Fig.  38  at  frequency  /  =  TOOTHz  .  (a) 
antisymmetric  mode,  (b)  symmetric  mode,  (c)  isolated  chain. 
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Figure  40  -  Amplitude  of  the  longitudinal  cross-polarization  for  the  chains  of  Figs.  36  and  38,  operating  in  the 

quasi-transverse  regime. 
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Figure  41  -  Analogous  to  Fig.  31  and  Fig.  36,  modal  dispersion  for  the  quasi-  transverse  z- polarized  modes 
supported  by  two  parallel  chains  with  interchain  distance  I  =  50  nm  . 
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Figure  13  -  (Color  online).  Similar  to  Fig.  12,  modal  dispersion  for  the  quasi-  transverse  z- polarized  modes 
supported  by  two  parallel  chains  with  interchain  distance  /  =  30  nm  . 
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Figure  41  shows  the  dispersion  of  symmetric  and  antisymmetric  transverse  modes  for  /  =  50  nm  . 
Here,  for  the  same  distance  as  in  Figs.  31  and  36,  the  coupling  perturbs  the  propagation 
properties  even  less  as  compared  to  the  isolated  chains. 

Also  in  this  case,  symmetric  modes  allow  slightly  longer  propagation  lengths  near  the  light  line, 
where  the  coupling  is  stronger.  Increasing  the  coupling  (/  =  30nm),  as  in  Fig.  42,  the 
perturbation  is  stronger,  even  if  the  trend  is  similar  as  in  the  previous  scenario.  Similar  to  the 
quasi-transverse  propagation  considered  in  the  previous  section,  these  purely  transverse  modes 
are  also  backward  in  nature  and  support  leaky-wave  propagation  in  the  upper  frequency  regime 
and  Bragg  stop-band  in  the  lower  one.  Moreover  in  this  case  the  anti- symmetric  transverse  mode 
(consistent  with  the  definition  in  Fig.  37)  supports  two  distinct  branches,  one  of  which  with  very 
weakly  guided  properties  near  the  light  line. 

Figure  43  shows  the  calculated  orthogonal  electric  field  distribution  in  the  xy  plane  for  the 
modes  of  Fig.  41  at  the  frequency  /  =  750  THz  .  In  this  case  the  modes  are  purely  transversely 
polarized  and  the  guided  wave  numbers  are  respectively  J3  =1.987-0.42/, 
Psym  =2.64-0.29 6/,  /FIlg|e  =2.36-0.34/,  consistent  with  Fig.  42.  The  field  confinement  in  this 

polarization  is  not  drastically  different  from  that  of  an  isolated  chain,  as  evident  from  the  figure, 
and  the  main  advantage  of  using  parallel  chains  may  reside  in  the  longer  propagation  distance  of 
symmetric  modes  near  the  light  line. 
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Figure  43  -  (Color  online).  Electric  field  distribution  (snapshot  in  time)  for  the  chains  of  Fig.  41  at  frequency 
/  =  750  77 /z  .  (a)  antisymmetric  mode,  (b)  symmetric  mode,  (c)  isolated  chain. 

e.  Conclusions 

We  have  presented  in  this  section  a  fully  general  and  complete  theoretical  formulation  for  the 
analysis  of  the  dynamic  coupling  between  two  parallel  linear  chains  of  plasmonic  nanoparticles 
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operating  as  optical  waveguides.  These  chains  may  support  up  to  eight  different  guided  modes 
with  different  polarization  properties  in  the  same  range  of  frequencies,  which  we  have  fully 
analyzed  here.  We  have  shown  that,  compared  to  linear  arrays,  these  waveguides  may  support 
relatively  longer  propagation  lengths  and  ultra-confined  beams,  operating  analogously  to 
transmission-line  segments  at  lower  frequencies.  In  particular,  our  results  confirm  that,  by 
operating  near  the  light  line  with  antisymmetric  quasi-longitudinal  modes,  we  may  achieve 
relatively  long  propagation  lengths  (of  several  wavelengths)  and  ultraconfined  beam  traveling, 
similar  to  a  transmission-line,  in  the  background  region  sandwiched  between  the  two 
antisymmetric  current  flows  guided  by  the  chains.  It  should  be  stressed  that  the  designs 
considered  here  are  based  on  ultrasmall  nanospheres,  with  the  aim  of  large  concentration  of  light 
in  a  sub-wavelength  region.  This  choice  inherently  produces  relatively  short  propagation  lengths, 
in  part  improved  by  the  parallel  chain  configuration.  Large  field  confinement,  even  at  the 
expenses  of  moderate  propagation  distances,  may  be  appealing  in  a  variety  of  applications,  i.e., 
nonlinearity  enhancement,  sensing,  optical  switching  and  nanoscale  interaction  with  light,  as  in 
optical  nanocircuits  [144]-[145].  On  the  other  hand,  longer  propagation  distances  may  be 
achieved  by  considering  larger  particles  or  lower  frequencies  of  operation  for  which  metals  are 
more  conductive,  as  discussed  for  isolated  chains  in  [140]  and  in  the  previous  sections.  We 
believe  that  these  results  may  be  of  great  interest  for  optical  connections  with  sub-wavelength 
lateral  cross-section  in  a  variety  of  new  generation  optical  devices  of  interest  to  U.S.  Air  Force. 


7.  Rigorous  First-Principle  Homogenization  Theory  of  Metamaterials 

a.  Summary 

We  introduce  in  this  section  a  first-principle  homogenization  theory  for  periodic  metamaterials 
composed  of  arbitrary  dielectric,  magnetic  and/or  conducting  inclusions.  We  derive  closed-form 
analytical  expressions  for  the  effective  constitutive  parameters  of  the  metamaterial  array, 
pointing  out  the  relevance  of  inherent  spatial  dispersion  effects,  present  even  in  the  long 
wavelength  limit.  Our  results  clarify  the  limitations  of  quasi-static  homogenization  models  and 
the  proper  physical  meaning  of  effective  metamaterial  parameters.  In  particular,  we  outline  the 
physical  reasons  behind  the  necessity  of  considering  magneto-electric  coupling  in  the 
homogenized  constitutive  relations,  even  in  the  case  of  metamaterials  formed  by  center- 
symmetric  inclusions.  These  results  may  be  of  great  interest  for  U.S.  Air  Force  applications  of 
nonconventional  materials  and  artificial  materials  in  realistic  systems  and  devices. 

b.  Introduction 

The  electromagnetic  homogenization  of  natural  and  artificial  materials  has  a  long-standing 
tradition  [  1 55]-[  1 60]  and  several  theories  are  available  to  properly  define  macroscopic  averaged 
quantities  representing  the  effective  constitutive  parameters  of  periodic  or  random  collections  of 
molecules  or  inclusions.  The  same  way  in  which  we  define  pennittivity  and  permeability  of 
natural  materials,  by  averaging  out  irrelevant  microscopic  fluctuations  of  the  fields  at  the  atomic 
or  molecular  level,  also  in  the  field  of  artificial  materials  and  mixtures,  homogenization  and 
mixing  rules  have  been  put  forward  over  the  years  in  order  to  avoid  solving  for  the  complex 
electromagnetic  interaction  among  a  large  number  of  inclusions  [160].  Homogenized 
descriptions  of  natural  and  artificial  materials  can  hold  only  in  the  long-wavelength  regime,  i.e., 
for  effective  wavelengths  and  averaged  field  variations  much  larger  than  the  material  granularity. 
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Within  these  limits,  however,  such  descriptions  have  proven  to  be  accurate  and  of  great 
advantage  for  analysis  and  design  purposes. 

With  the  advent  of  metamaterials  [161],  i.e.,  artificial  materials  with  anomalous  and  exotic 
electromagnetic  response,  the  necessity  of  more  advanced  concepts  and  improved 
homogenization  models  has  become  evident,  since  often  the  topological  and/or  resonant 
properties  of  the  inclusions  and  building  blocks  do  not  allow  a  description  of  their  properties  in 
terms  of  simple  averaging  procedures.  Often,  the  exotic  metamaterial  properties  are  inherently 
based  on  these  anomalous  features,  which  cannot  be  captured  by  simple  homogenization 
schemes  inspired  to  natural  materials.  The  necessity  for  improved  homogenization  models  has 
been  outlined  in  several  recent  papers  on  the  topic  [162]-[179],  which  describe  different 
approaches  to  the  problem. 

The  simplest  homogenization  method  is  represented  by  retrieval  techniques,  which  postulate  the 
equivalence  between  a  complex  metamaterial  array  and  a  unifonn  slab  of  same  thickness  with 
unknown  constitutive  parameters,  often  limited  to  pennittivity  and  permeability  [  1 80]-[  1 82] . 
This  simple  approach  regularly  provides  non-physical  dispersion  of  these  constitutive 
parameters,  in  particular  near  the  inclusion  resonances,  yielding  complex  values  of  permittivity 
and  permeability  that  violate  basic  passivity  and  causality  constraints  [183].  The  typical  presence 
of  “anti-resonant”  artifacts  in  the  dispersion  of  the  effective  constitutive  parameters,  wrong  sign 
of  the  imaginary  part  and  wrong  slope  of  the  real  part  of  the  extracted  parameters  are  all  clear 
signs  of  the  inadequacy  of  such  simple  homogenization  approaches  when  applied  to  resonant 
artificial  materials,  as  discussed  in  [184], 

These  anomalies  have  been  generically  related  to  strong  spatial  dispersion  effects  present  in 
metamaterials,  which  should  be  taken  into  account  in  more  refined  homogenization  models.  In 
this  context,  analytical  and  semi-analytical  methods  have  been  put  forward  to  address  the 
homogenization  in  a  more  rigorous  fashion.  Generalized  Clausius-Mossotti  techniques  have  been 
extended  to  the  case  of  complex  inclusions,  including  bianisotropic  effects,  possible  presence  of 
spatial  dispersion  and  accurate  modeling  of  the  interaction  among  inclusions  [  1 67]-[  175],  As 
another  successful  approach,  averaging  a  planar  sheet  of  inclusions  and  then  considering  the 
mutual  interaction  among  parallel  layers  as  a  Bloch  lattice  has  been  proposed  as  a  venue  to 
derive  physically  sound  effective  constitutive  parameters  [  162]-[1 66].  In  these  schemes  too, 
however,  effects  of  spatial  dispersion  can  often  generate  artifacts  in  the  extracted  effective 
parameters,  which  are  not  easily  explained  on  physical  grounds. 

In  order  to  circumvent  these  issues,  an  interesting  approach  to  the  homogenization  of  periodic 
arrays  of  dielectric  inclusions  [176]- [179]  has  been  put  forward  based  on  the  model  that 
commonly  describes  optical  crystals  [185],  i.e.,  introducing  a  generalized  permittivity  tensor  that 
includes  all  the  polarization  effects,  including  artificial  magnetism  and  possible  bianisotropy. 
Although  this  technique  is  applicable  to  a  wide  class  of  dielectric  metamaterials  with  periodic 
features,  it  often  may  not  be  intuitive  to  relate  weak  spatial  dispersion  effects  in  the  generalized 
permittivity  tensor  to  magnetic  or  bianisotropic  effects,  and  it  may  be  more  desirable  to  describe 
these  effects  in  terms  of  local  permeability  or  chirality  parameters.  In  addition,  in  this  description 
the  generalized  permittivity  tensor  explicitly  depends  on  the  wave  number  even  for  local 
materials,  making  it  more  challenging  to  apply  the  usual  boundary  conditions  and  the  solution  of 
dispersion  relations. 

In  most  circumstances,  it  may  be  safely  assumed  that  for  frequencies  well  below  the  inclusion 
resonances,  for  which  both  the  background  wavelength  and  the  effective  guided  wavelength  are 
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significantly  larger  than  the  average  array  granularity,  local  or  quasi-local  constitutive 
parameters  should  be  sufficient  to  characterize  the  metamaterial  response,  and  strong  spatial 
dispersion  effects  may  be  negligible.  There  are,  however,  examples  of  metamaterial  geometries 
that  contradict  this  general  assumption,  such  as  the  wire  medium  [186]  or  dense  arrays  of 
plasmonic  nanoparticles  [169],  [176],  which  show  spatial  dispersion  effects  even  for  longer 
wavelengths.  In  addition  to  these  effects,  anomalous  bianisotropic  coupling  due  to  lattice  effects 
has  been  pointed  out  in  recent  papers  [174]-[175]. 

Inspired  by  all  these  issues,  in  the  following  we  develop  a  general  and  rigorous  first-principle 
theory  for  the  homogenization  of  periodic  metamaterials  composed  of  arbitrary 
magnetodielectric  and/or  conducting  inclusions.  We  use  a  homogenization  approach  consistent 
with  [176], [185],  but  generalized  to  the  case  of  an  array  comprised  of  arbitrary  electric  and 
magnetic  materials  and  to  the  presence  of  arbitrary  sources.  We  show  that  it  is  indeed  possible  to 
derive  a  fully  consistent  homogenization  model  for  metamaterial  arrays,  which  rigorously  takes 
into  account  their  complex  wave  interaction,  that  does  not  depend  on  the  external  source 
excitation  and  that  converges  to  quasi-local  constitutive  parameters  in  the  long-wavelength  limit. 
Our  analysis  highlights  the  reasons  and  physical  mechanisms  behind  artifacts  and  non-physical 
features  arising  in  common  homogenization  models,  showing  that  a  rigorous  analysis  of  the  array 
coupling  inherently  requires  considering  frequency  and  spatial  dispersion  effects  at  the  lattice 
level,  even  when  higher-order  multipolar  interaction  or  bianisotropic  effects  within  the  unit  cell 
are  negligible.  These  effects  modify  the  usual  fonn  of  effective  local  constitutive  parameters  and 
they  may  be  associated  with  a  direct  manifestation  of  the  finite  phase  velocity  with  respect  to  the 
array  granularity,  becoming  particularly  relevant  for  more  densely  packed,  and  possibly  resonant, 
metamaterials.  Our  findings  establish  the  foundations  of  a  physically  meaningful  and  rigorous 
description  of  a  wide  class  of  metamaterials,  in  particular  when  the  density  of  inclusions  is  not 
small  and  classic  homogenization  models,  like  Clausius-Mossotti  relations  [160],  lose  their 
accuracy  even  in  the  long  wavelength  regime.  In  addition,  we  shed  new  light  into  the  physics  of 
metamaterials  and  the  meaning  of  effective  constitutive  relations  to  describe  their  complex  wave 
interaction. 


c.  General  Homogenization  Theory  for  Metamaterial  Arrays 

In  this  section,  starting  from  first-principle  considerations,  we  develop  a  general  and  rigorous 
homogenization  model  for  periodic  arrays  of  arbitrarily  shaped  dielectric,  magnetic  and/or 
conducting  inclusions,  extending  the  approach  commonly  used  in  optical  crystals  [185]  and 
dielectric  metamaterials  [176]  to  arbitrary  materials  and  arbitrary  fonn  of  excitation.  We  will 
show  that  this  general  description  reduces  to  a  quasi-local  constitutive  representation  in  the  long- 
wavelength  limit.  For  simplicity  of  notation,  we  assume  a  cubic  lattice  with  period  d ,  but 
extension  to  arbitrary  lattices  may  also  be  envisioned.  The  most  general  description  of  a  periodic 
array  in  its  linear  operation  may  be  developed,  without  loss  of  generality,  in  the  Fourier  domain 
[185].  Our  goal  is  to  derive  the  general  fonn  of  macroscopic  constitutive  relations  for  any 
arbitrary  pair  ( to,  [I ) ,  relating  spatially  averaged  field  quantities  that  vary  as  e^  re~,<ot .  In  general, 
only  a  limited  set  of  eigen-vectors  p  are  supported  by  the  anay  at  a  given  frequency  to  in 
absence  of  impressed  sources.  These  correspond  to  the  eigen-modes  of  the  system.  Therefore,  we 
will  assume  the  presence  of  impressed  sources  with  specific  e^'re~,wt  plane -wave  like 
dependence,  uniformly  distributed  all  over  the  anay.  This  ensures  an  averaged  space-time 
distribution  of  the  induced  fields  with  the  same  e'pre"“"  dependence.  In  practice,  it  may  be 
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challenging  to  practically  realize  a  distribution  of  uniformly  impressed  sources  within  a 
metamaterial  array,  so  this  excitation  should  be  seen  as  a  test  excitation  to  isolate  the  response  in 
the  Fourier  domain,  or  as  the  Fourier  expansion  of  embedded  sources  with  realistic  space-time 
distribution.  We  will  specialize  these  results  to  the  eigen-modal  scenario  (source-free  modes) 
later  in  the  paper. 

In  the  most  general  case,  the  microscopic  field  distribution  [187]  at  any  point  in  the  array 
satisfies 

VxE(r)=»ftH(r)ti»M(r)- K„e»' 

,  (46) 

v  x  H  (r)  =  -i(O£0 E  (r )  -  id. P  (r )  +  Jexte‘l>  r 

where  E(r),  H(r)  are  the  local  electric  and  magnetic  fields,  P(r)  is  the  local  polarization 
vector,  M(r)  is  the  local  magnetization  vector,  Jext  and  Kext  are  complex  vectors  of 
independently  impressed  distributed  electric  and  magnetic  current  density  sources  with  explicit 
plane  wave  dependence  e,f'T ,  and  s0 ,  ju0  are  the  background  permittivity  and  penneability, 

repectively.  Due  to  the  linearity  of  the  problem,  we  have  suppressed  in  (46)  the  e~l<at  time 
dependence.  In  the  presence  of  electric  or  magnetic  conductors,  the  induced  current  densities 
Jind ,  K.nd  are  implicitly  embedded  into  P (r)  =  iJind  (r) /  co  and  M(r)  =  z'KiW  (r)  /  0)  in  Eq.  (46) 


The  distributed  impressed  source  distributions  may  also  be  seen  as  sustaining  impressed  fields 
with  the  same  e‘p re~uot  plane-wave  dependence  and  therefore  complex  amplitudes  satisfying 

*P  X  -  KeU  (4y) 

/PxHm=-/^oEe„+J„/ 


Notice  that  the  arbitrary  choice  of  Jew  and  K  ext  in  (46)  implies  that  the  complex  amplitudes  of 
impressed  fields  Eex(  and  H  w  are  independent  of  each  other.  This  will  be  very  important  to 

ensure  the  general  validity  of  the  effective  homogenization  model  proposed  here,  as  discussed 
below. 


Due  to  the  periodicity  of  the  crystal,  we  may  write  Eq.  (46)  in  the  e'pre  Fourier  domain 


z'P  x  E  =  icoju0 H  +  icoM  -  K(,v; 
/(I  x  H  =  -icos0  E  -  icoP  +  Jext 


(48) 


where  the  bar  denotes  the  averaging  operation  E 


— -J( E(r)e  ,prdr,  and  similarly  for  all  the 


other  vectors  in  Eq.  (48).  This  averaging  procedure,  consistent  with  [185],  [176],  filters  out  the 
dominant  contribution  to  the  local  field  E(r) ,  varying  as  Ee'|lr ,  of  interest  for  a  macroscopic 


homogenized  description  of  the  array.  Eq.  (48)  effectively  relates  the  complex  amplitudes  of  the 
spatially  averaged  macroscopic  field  quantities,  which  all  vary  with  an  e'pre~wr  space-time 
dependence  due  to  the  chosen  form  of  impressed  excitation  and  the  linearity  of  the  problem. 
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Inspecting  Eq.  (48),  one  may  be  tempted  to  define  spatially  averaged  displacement  vectors 
B  =  //0H  +  M  and  D  =  £0E  +  P  ,  and  the  associated  constitutive  relations  B  =  •  H ,  D=£^-E, 

which  would  generalize  the  metamaterial  homogenization  approach  used  in  [185],  [176]  to  the 
case  of  magnetodielectric  materials.  However,  this  macroscopic  description  would  have  several 
shortcomings:  the  effective  pennittivity  and  penneability  coincide  with  s0  or  //0  if  the  inclusions 

are  formed  by  purely  magnetic  (P(r)  =  P  =  0)  or  dielectric  (M(r)  =  M  =  0)  materials, 

respectively.  This  implies  that  artificial  magnetic  or  polarization  effects,  stemming  from  the 
rotation  of  electric  and  magnetic  polarization  respectively,  remain  hidden  as  spatial  dispersion 
effects  in  the  pennittivity  or  penneability  tensors.  In  particular,  coincides  with  the 

generalized  permittivity  defined  in  [176]  in  the  case  of  dielectric  inclusions.  This  description, 
therefore,  cannot  converge  to  a  local  constitutive  model  in  the  long-wavelength  limit  in  the 
presence  of  artificial  magnetic  or  dielectric  effects  (e.g.,  for  split-ring  resonators),  of  most 
interest  for  metamaterial  applications. 


i.  Multipolar  expansion 

In  order  to  overcome  this  problem,  we  assume  that  the  unit  cell  is  sufficiently  smaller  than  the 
wavelength  of  operation  to  ensure  that  the  induced  microscopic  polarization  and  magnetization 
vectors  vary  slowly  within  each  unit  cell.  In  such  circumstances,  it  is  possible  to  expand  P  in 
Taylor  series  around  the  origin  of  each  unit  cell,  to  obtain  [188] 


-If  P(r) 


e'^dr  = 


1 


J,P(r)^r  +  ?'Pxj/X^r)^-f-J,[rP(r)  +  P(r)r]jr 
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IV  2 

rxP(r)xr 
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dr  x  p  +  ^  x 


H 


[rxP(r)]r  +  r  [rxP(r)] 
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+  . 


pxM£  7p  ^  ,  P 
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P(r) 
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represent  the  first  electric  and  magnetic  multipole  moments  associated  with  the  induced  electric 
polarization  distribution  P(r).  In  particular,  P;  ,  M£  are  the  fist-order  contribution  to  the 

electric  and  magnetic  dipole  moments,  respectively,  QE ,  Q"'  are  the  electric  and  magnetic 
quadrupole  moment  cntributions,  P)  is  the  third-order  contribution  to  the  electric  dipole 
moment.  The  subscript  E  indicates  the  origin  of  these  multipole  moments,  associated  to  the 
electric  polarization.  We  can  apply  analogous  considerations  to  the  induced  magnetization  M(r) 


M  =  Mh  + 


PxP„  /p  ^  _  p 


0)£n 


f-Q  H+->< 
2  ~  ico 


x  p  + 


pq; 


+... , 


with  analogous  definitions  for  the  corresponding  multipole  moments 

M*=-UM(r)rfr 


d' 

ico£0 


I, 


r  x 


M(r) 


dr 


Q  H  =  \v  [rM  (r)  +  M  (r  )  r]  dr 
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,  _  id)  r 


io)  r  r xM(r)xr 


dr 
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ICO 
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[rxM(r)]r  +  r[r  xM(r)] 


dr 


v  3 

and  the  subscript  H  refers  to  the  magnetic  origin  of  these  multipole  moments. 
Using  the  previous  expansions,  we  may  write  Eq.  (48)  in  the  interesting  form 


zpx 

Zpx 


E_Z»_,pxM'„  +  i.Q; 

£n  2 


H 


'0 

M 


=  10) 


//0H  +  M,-|.Q“ 


/P 


-ZPxP;-^ 

Mo  2  - 
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=  -ico 

^0e+p£-^.q; 

x 
J 
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K 


(51) 


(52) 


(53) 


in  which  we  have  neglected  the  effects  of  higher  order  multipole  moments  (beyond  the  electric 
and  magnetic  quadrupole  moments)  in  Eqs.  (49),  (51). 


ii.  Proper  definition  of  averaged  fields 

Eq.  (53)  ensures  that,  by  correcting  the  definition  of  the  averaged  fields  as  follows: 
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E,„.=E-i-,pxM>iQ‘„ 

^0  2  ~ 

h„-h-m' 

A,  2  _ 

? 

D„=^+rf“Q'f 

the  macroscopic  (averaged)  Maxwell’s  equations  take  the  expected  fonn 
/PxEav  =/ABav-K^ 

/PXHav  =-^D,v+Je« 


(54) 


(55) 


for  any  arbitrary  pair  (ry,p) . 

The  definition  of  averaged  fields  as  given  in  (54)  solves  the  issues  outlined  above  in  reference  to 
Eq.  (48)  and  ensures  the  proper  representation  of  artificial  electric  and  magnetic  effects,  making 
sure  that  the  constitutive  relations  of  the  metamaterial  tend  to  local  parameters  in  the  long- 
wavelength  limit,  even  in  presence  of  artificial  magnetic  or  polarization  effects.  In  other  words, 
this  averaging  procedure,  based  on  a  rigorous  first-principle  approach,  properly  takes  into 
account  weak  forms  of  spatial  dispersion  associated  with  artificial  magnetism  and  polarization,  at 
the  basis  of  common  metamaterial  effects,  and  it  allows  their  description  in  a  local  framework. 
Eq.  (54)  shows  that  the  proper  form  of  averaged  electric  and  magnetic  fields  EflV  and  Hm,  is 

obtained  after  correcting  the  spatial  microscopic  averages  E  and  H  for  the  possible  presence  of 
these  artificial  effects,  associated  with  the  rotation  of  P  and  M .  This  ensures  that  they  are 
correctly  attributed  to  local  constitutive  parameters  in  the  long-wavelength  limit.  In  the  presence 
of  only  dielectric  and  conducting  inclusions  (M(r)  =  0),EflV=E  and  BflV  =  //0H  ,  ensuring  that 

E  and  B  are  the  direct  source  fields,  consistent  with  rigorous  homogenization  for  optical 
crystals  and  dielectric  metamaterials  [185], [176],  Similarly,  if  only  magnetization  is  present  at 
the  microscopic  level,  then  H  av  =  H  and  Da  =  s0 E  are  the  source  fields  to  be  spatially 
averaged,  as  shown  in  [178],  Eq.  (54)  represents  a  generalization  of  these  techniques  in  the  case 
of  magnetodielectric  arrays,  for  which  both  averaged  fields  Eav ,  Ha  need  to  be  corrected  for  the 

possible  presence  of  artificial  electric  and  magnetic  effects.  This  is  the  only  way  to  ensure  that 
these  weak  spatial  dispersion  effects  are  properly  taken  into  account  within  a  local 
homogenization  description. 


Hi.  Relations  between  averaged  fields 

In  the  general  case,  the  constitutive  relations  relating  the  averaged  displacement  vectors  Dav , 
Ba  to  the  averaged  field  vectors  Eav,  Hav  inherently  depend  on  p,  implying  that  spatial 
dispersion  effects  associated  with  higher-order  multipole  contributions  are  present.  In  the  long 
wavelength  limit  of  interest  here,  however,  it  is  expected  that  the  distributions  of  P(r) ,  M(r) 
may  be  described  simply  in  terms  of  electric  and  magnetic  dipole  moments,  which  is  the  case 
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when  the  unit  cell  is  sufficiently  smaller  than  the  wavelength  of  operation,  and  the  inclusions  are 
not  too  densely  packed.  In  such  circumstances,  the  explicit  effects  of  spatial  dispersion  disappear 
in  Eq.  (54): 

E«v  =E-Pff /fi0 

Hav=H-M£///0 

(56) 

Dflv  -  *0E  + 

®av  = 

and  the  constitutive  relations  relating  averaged  displacement  and  field  vectors  may  be  written  in 
the  local  fonn 


D  =  £•  E  +  P„  +  P„  =  en  E  +  P 

av  0  av  E  H  0  av  av 

B  =  «„H  +  M„+M  =«„H  +M 

av  r'O  av  H  E  r*0  av  av 


(57) 


where  we  have  combined  averaged  polarization  and  magnetization  stemming  from  electric  and 
magnetic  effects  into  Pav  and  M  av .  As  an  aside,  this  relation  rigorously  proves  from  first- 

principle  considerations  that  the  electric  and  magnetic  contributions  are  simply  additive  within 
the  unit  cell  in  the  dipolar  limit,  which  in  turn  implies  that  we  can  alternatively  describe  the 
complex  unit  cell  interaction  in  terms  of  point-dipoles.  In  other  words,  the  following  theory  may 
be  also  obtained  in  a  more  streamlined  version  by  considering  lattices  of  electric  and  magnetic 
point-dipoles  to  describe  the  complex  inclusions  in  each  unit  cell. 

Combining  Eq.  (57)  and  (47)  into  (55),  we  can  write  a  general  relation  between  averaged  and 
external  fields  and  averaged  polarization  and  magnetization  vectors, 

/p  x (Em,  - Eext )  =  toft,  (H„  -  H  , )  +  icoMav 

'P  X  (Hav  -  H  )  =  -i®£0  (  E«,  -  Eexr  )  -  ^Pav 

These  equations  may  be  further  manipulated  to  yield 

[kl  +  p  x  p  x]  (Eav  -  Eext )  =  -k;  ^  +  k0770P  x  ^ 

Ac 

[C02  +  P  x  p  x] (Hav  -  )  =  -kl  -  —  P  x  — 

Ao  % 

where  k0  =  a>^£0ju0  and  //0  =  ^ju0  /  £0  . 


(58) 


(59) 


Henceforth,  for  simplicity  of  notation,  we  consider  only  averaged  and  impressed  field 
distributions  that  are  transverse-electromagnetic  (TEM)  waves  propagating  along  the  direction  p 
(where  the  hat  indicates  a  unit  vector).  It  may  be  proven  that  this  is  the  case  for  isotropic  unit 
cells  and  lattices,  as  we  assume  in  the  following,  when  p  is  along  one  of  the  lattice  axes  or,  more 
generally,  for  any  propagation  direction  in  the  limit  k0d  «:  1 .  A  more  general  tensorial  notation 

may  be  used  for  arbitrary  propagation,  but  is  not  adopted  here  in  the  interest  of  notational 
simplicity.  In  this  situation,  Eq.  (59)  may  be  compactly  written  as  follows, 
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This  is  a  very  general  result,  which  relates  the  averaged  and  external  fields  for  any  (&>,p)  pair 

and  holds  for  any  metamaterial  array  and  any  combination  of  electric  and  magnetic  excitations. 
Observe  that,  analogous  to  the  way  both  electric  and  magnetic  induced  currents  contribute  to  the 
electric  and  magnetic  averaged  polarization  effects  in  (57),  both  averaged  polarization  and 
magnetization  vectors  contribute  to  the  averaged  electric  and  magnetic  fields.  In  other  words,  an 
inherent  form  of  magnetoelectric  coupling  at  the  unit  cell  level  stems  from  weak  spatial 
dispersion  effects  when  J3  *  0 ,  associated  with  finite  phase  velocity  across  each  unit  cell.  Eq. 
(60)  defines  a  general  relation  among  averaged  and  impressed  fields,  which  is  independent  on 
the  specific  nature  of  the  metamaterial  inclusions.  In  the  following  section,  we  use  these  results 
and  introduce  the  inclusion  into  the  picture,  in  order  to  derive  the  first-principle  effective 
constitutive  parameters  of  an  arbitrary  metamaterial. 


iv.  Effective  constitutive  parameters 

After  having  established  the  proper  definition  of  averaged  fields  and  their  general  relations,  we 
are  now  ready  to  derive  a  macroscopic  homogenized  description  of  the  array,  once  we  relate 
averaged  polarization  and  magnetization  vectors  to  the  local  fields,  as  a  function  of  the  specific 
inclusion  geometry.  Since  we  are  assuming  that  dipolar  terms  are  dominant,  we  may  compactly 
describe  the  unit  cell  response  in  terms  of  its  polarizability  coefficients,  which  relate  the  induced 
electric  and  magnetic  dipole  moments  in  the  unit  cell  p000  =  d?'Pnv  and  m000  =  d"'M av  to  the  local 
fields  at  its  center. 


Pooo  =  £0aeE,oc  -  X  H/oc 


m000  =MoamRloc-Moae, 


PxE 


loc 


(61) 


Vo 


where  ae ,  am  and  aem  are  the  electric,  magnetic  and  magneto-electric  polarizability  coefficients, 

respectively.  All  these  coefficients  have  dimensions  of  an  inverse  volume,  and  they  are 
considered  scalar  here  due  to  the  assumptions  of  TEM  propagation  and  isotropic  unit  cell.  In 
addition,  in  writing  Eq.  (61)  we  have  implicitly  assumed  that  the  inclusions  are  reciprocal. 

The  fields  E/oc  and  H/oc  represent  the  local  fields  at  the  center  of  the  unit  cell  in  absence  of  the 
inclusion.  They  are  due  to  the  superposition  of  the  impressed  fields  Eext ,  H  ,a  and  the  induced 
fields  scattered  from  inclusions  in  the  rest  of  the  array,  that  is, 
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The  interaction  constants  C  and  Cem  may  be  evaluated  using  the  dipolar  radiation  from  the 
generic  unit  cell  at  xlmn-{ld,md,nd )  and  applying  the  Floquet  condition  plmn  =  p000elfi'rimn  , 

m  =  m  p'$'rlmn  • 

111  Imn  111 000 1 ^ 

C=  X  P-Q,Ar,mn)^' .  P 

{l,m,m)*(  0,0,0) 

Q»  =  X  P'Gem  (»*/,„„  )e'pr"”  m 

(/,m,m)^(0,0,0) 


where  Gee(r/mJI)  and  Gem  (r/mn)  are  the  electric  and  magneto-electric  dyadic  Green’s  functions 
[189]  and  p,  m  are  unit  vectors  oriented  along  p000  and  m000.  Fast  converging  expressions  for 
these  summations  are  available  in  [168]-[169],  [173], 

Combining  Eqs.  (6 1  )-(62),  we  may  now  derive  a  general  relation  between  impressed  fields  and 
averaged  polarization  vectors: 
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which,  substituted  in  (60),  provides  the  important  relations 


E  = 


H 


d\x„ 


CX  CX  4-  CX 


2  ~dCim 


-  + 


d  a, , 


a  a  +  a 

e  m  em 


2 


°o 

M 


<73a„ 


✓y  /y  -f  CX 


-  +  J3C 


3^ 
em 


%P: 


M 


Ao 


-  + 


Ao 


d  a „ 


a  a  +  a 

e  m  em 


em 


- d3C ; 


P  P 

—  X  — 

7o  *0 


(65) 


Here  we  have  used  the  reduced  interaction  constants 


C,=C- 


C'  =C 

em  em 
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d  p2-k l 

i  PK 

d 3  p2-k 2 


(66) 


which  respectively  coincide  with  p  •  Cm,  p  and  p  C  m  m  derived  in  [176]  using  an  alternative 
spectral-domain  representation. 

Eq.  (65)  represents  another  important  result,  since  it  shows,  directly  from  first-principle 
considerations,  that  it  is  possible  to  establish  a  general  relation  between  averaged  electric  and 
magnetic  fields  and  averaged  polarization  vectors  [as  defined  in  (56)],  which  depends  on  the 
array  period  and  the  polarizability  coefficients  for  any  given  pair  (&>,  p) ,  but  is  completely 
independent  on  the  relative  amplitude  of  the  impressed  sources  Jext ,  Kext .  This  is  a  fundamental 
property  of  this  homogenization  theory,  since  a  proper  description  of  the  metamaterial 
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interaction  in  homogenized  form  cannot  depend  on  the  type  and  form  of  external  excitation,  as 
commonly  happens  in  less  rigorous  homogenization  models. 

The  relations  in  (65)  also  show  that  there  is  an  inherent  form  of  magnetoelectric  coupling 
(usually  negligible  in  natural  materials)  relating  Em,  to  the  rotation  of  M  av  and  Hav  to  the 

rotation  of  Pav .  As  expected,  part  of  this  coupling  is  associated  with  the  presence  of  aem ,  which 

represents  the  possible  bianisotropy  within  the  unit  cell,  stemming  from  asymmetric  or  non- 
centered  inclusions  [191].  However,  Eq.  (65)  predicts  that,  even  when  inclusions  are  perfectly 
center-symmetric  and  with  no  inherent  bainisotropy,  a  form  of  magneto-electric  coupling  is  still 
expected,  associated  with  the  presence  of  C'em .  This  additional  coupling  term  is  due  to  lattice 
effects  and  the  finite  value  of  [3  .  We  will  discuss  the  implications  of  this  coupling  in  more  detail 
in  the  following. 

Using  (57),  we  can  finally  write  for  the  constitutive  relations  of  the  metamaterial  array: 

D  =  £nE  +  P  =  s  „E  -  ( ye„.  +  y°„ )  B  x  H 

av  0  av  av  eff  av  \  A  eff  XI  eff  J  r  av 

,  ^  >  (67) 

v  =  Ao  Hav  +  Mav  =  pieyWav  —yZeff  —  Xeff  )  P  X  ^av 
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where  c0  =  1  /  4^ ,  a4ff)  =  a>" 1 A  >  a4ff)  =  a* 1 A  ’  a^W)  =  a,n  /  A  and  A  =  (aeam  +  a2em )  [in 

the  absence  of  magnetoelectric  coupling  at  the  unit  cell  level  aem  =  0  and  ae(cjf)  =  ae, 
a  ,  =  a„  1. 

myeff)  m  J 


v.  General  properties  of  the  effective  constitutive  parameters 
These  expressions  represent  general  closed-fonn  effective  constitutive  parameters  obtained  from 
first-principle  considerations.  They  are  valid  for  any  pair  (ru,p)  and  any  form  of  external 

excitation  Jext ,  Kext ,  ensuring  that  this  homogenized  description  does  not  depend  on  the  specific 

impressed  field  distribution  in  each  unit  cell,  but  it  is  the  inherent  response  of  the  metamaterial  as 
a  bulk  to  a  given  arbitrary  level  of  electric  and/or  magnetic  excitation.  It  is  noticed,  in  particular, 
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that  the  ratio  of  averaged  fields  Eav  /  Hav ,  i.e.,  the  local  wave  impedance,  inherently  depends  on 
the  specific  choice  of  impressed  sources  Jexl ,  Kext ,  following  Eq.  (58),  as  expected  in  presence 

of  arbitrary  impressed  sources.  However,  the  constitutive  parameters  defined  in  Eq.  (68)  do  not 
depend  on  this  ratio,  thus  compactly  including  the  macroscopic  polarization  and  magnetization 
properties  of  the  array  for  arbitrary  excitation. 

Effective  permittivity  and  permeability  are  found  in  closed  form  in  the  first  two  expressions  (68). 
These  generalize  the  Clausius-Mosotti  homogenization  formulas  [  1 59]-[  1 60],  [176]  by 
rigorously  taking  into  account  the  coupling  among  the  inclusions  and  their  polarization 
properties.  More  importantly,  this  theory  demonstrates  the  inherent  presence  of  the  magneto¬ 
electric  coupling  via  the  coefficients  xljf  ar|d  X°eff  in  Eq.  (67)-(68).  The  first  portion  of  the 

bianisotropy  coefficient  xljf ,  even  with  respect  to  [3 ,  is  associated  with  the  magneto-electric 

effects  at  the  inclusion  level,  and  satisfies  the  usual  reciprocity  constraints  for  bulk  materials.  An 
additional  contribution  to  xe/f  >  °dd  with  respect  to  fi  ,  is  associated  to  inherent  magneto-electric 

coupling  effects  at  the  lattice  level.  These  latter  effects  cannot  be  generally  neglected.  Even  in 
the  case  of  center-symmetric  inclusions,  for  which  aem  =  0 ,  this  latter  contribution  is  present  as 

long  as  C'em  ^  0 .  The  presence  of  this  odd  bianisotropy  coefficient  has  been  pointed  out 
theoretically  and  numerically  in  [  1 74]-[  175],  and  the  present  theory  explains  its  physical  nature 
and  relevance  from  first-principle  considerations:  x°ejf  is  a  weak  spatial  dispersion  effect 

associated  with  the  finite  phase  velocity  of  mode  propagation  along  the  array,  not  negligible  even 
in  the  long-wavelength  limit  as  we  show  in  the  following  numerical  examples  (Section  8).  Its 
nature  is  associated  with  the  inherent  asymmetry  introduced  by  the  phase  propagation  across  the 
unit  cell,  and  this  explains  why,  at  first  sight,  its  occurrence  in  (67)  does  not  appear  to  satisfy  the 
reciprocity  relation  for  local  bianisotropic  materials.  Its  odd  response  with  respect  to  f3  ensures 
however  that,  by  reversing  the  direction  of  propagation  for  given  frequency,  its  contribution  also 
changes  sign,  ensuring  that  the  constitutive  relations  (67)  are  effectively  describing  a  reciprocal 
medium.  This  shows  the  drastic  difference  between  the  bianisotropy  stemming  from  lattice 
effects  x°eff  ar|d  the  one  associated  with  an  actual  magneto-electric  coupling  at  the  inclusion  level 

Xlff  ■  Its  relevance  in  the  homogenization  of  metamaterials  and  in  restoring  the  physical  meaning 
of  their  constitutive  parameters  has  been  discussed  in  further  detail  in  [184], 

Due  to  the  inherent  properties  of  the  summations  in  (63),  which  for  real  f3  satisfy  [168]-[169], 
[173],  [176] 


Im[C]  =  k\  /  (6;r) 
Im[C„,]  =  0 


(69) 


and  the  lossless  conditions  on  the  polarizability  coefficients  [190] 
hnfa;1]  =  lm[crj]  =  k\  I{6tt) 

|mKJ=° 

it  is  recognized  that  all  the  constitutive  parameters  in  (68)  are  real  for  lossless  particles  and  real 
1 3  ,  as  required  in  lossless  bianisotropic  materials,  satisfying  power  conservation. 
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Before  concluding  this  section,  it  is  worth  stressing  that  the  closed-form  expressions  (68)  for  the 
effective  parameters  of  the  array  apply  to  any  plane-wave  like  field  distribution  in  the 
homogenized  material,  any  form  of  excitation  and  any  pair  (&>,  p) ,  representing  an  accurate  first- 

principle  homogenization  model  for  the  metamaterial  array.  The  derived  parameters  are,  in  the 
general  case,  still  weakly  dependent  on  /? ,  as  a  symptom  of  spatial  dispersion.  However,  they 
tend  to  local  parameters  in  the  long- wavelength  limit  (small  /3 ).  The  present  general  theory 
highlights  that,  in  addition  to  artificial  magnetism  and  polarization  effects  stemming  from  weak 
spatial  dispersion,  the  rigorous  consideration  of  the  coupling  within  the  array  requires  additional 
magneto-electric  coefficients,  even  in  the  case  of  center-symmetric  inclusions.  In  the  following 
section,  we  consider  the  special  case  of  eigen-modal  solution,  without  external  impressed 
sources. 


d.  Eigen-Modal  Propagation  and  Equivalent  Constitutive  Parameters 
In  the  eigen-modal  case,  i.e.,  in  the  absence  of  external  distributed  sources,  Eq.  (64)  ensures  that 
non-trivial  solutions  are  available  only  for  specific  instances  of  p(&>),  satisfying  the  array 
dispersion  relation 

=  (71) 

The  corresponding  eigenvectors,  again  solving  Eq.  (64),  satisfy 

PoQO'P  __  P,  _  '  ^em+aem(eff)  _  1  a>n(eff)  ^ 

mooo-m  Mav  rj0  a^-C  Jj0  Cem  -  aem^eff)  ’ 

which  provides  a  specific  constraint  on  the  ratio  Pav  /  Mav  in  this  source-free  case,  function  only 
of  the  frequency  and  the  array  geometry. 

Rearranging  Eq.  (67)  and  (58),  in  this  regime  we  may  also  write 
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E  =  - icos  E 

dv  eq  av 


(73) 


which  shows  that  the  eigen-modal  propagation  may  be  described  in  terms  of  equivalent 
permittivity  and  penneability  parameters  seq ,  p  ,  which  embed  the  magneto-electric  coupling 

effects  as  a  form  of  weak  spatial  dispersion.  Their  validity  is  strictly  limited  to  eigen-modal 
propagation,  since  the  ratio  Pav  /  Mav  is  in  general  a  function  of  the  impressed  distributed 
sources.  The  description  of  the  array  in  terms  of  equivalent  parameters  is  particularly  attractive 
in  the  absence  of  bianisotropic  effects  at  the  inclusion  level  ( aem  -  %eeff  -  0 ),  for  which  the 

residual  magneto-electric  coupling  associated  with  lattice  effects  may  be  embedded  into 
equivalent  permittivity  and  permeability  related  to  the  effective  parameters  through  the 
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c  y 

normalization  factor  1 —  .  This  shows  that  classic  homogenization  models  that  aim  at 

j3/k0 

describing  metamaterial  arrays  in  tenns  of  pennittivity  and  penneability  (see,  e.g.,  [  1 66]-[  1 73]) 
effectively  extract  these  equivalent  parameters,  that  inherently  contain  a  fonn  of  weak  spatial 
dispersion  when  %°r  is  not  negligible.  It  is  evident  that  this  may  easily  translate  into 

inconsistencies  and  lack  of  physical  meaning  in  the  extracted  or  retrieved  parameters,  as 
discussed  more  in  detail  in  [184].  It  is  also  worth  stressing  that  these  parameters  are  inherently  a 
function  of  the  specific  ratio  Pav  /  Mav  in  (72),  i.e.,  they  are  forced  to  change  when  impressed 

sources  are  present  in  the  array,  in  sharp  contrast  with  the  effective  constitutive  relations,  derived 
in  (68)  on  first-principle  physical  grounds. 


i.  Secondary  parameters  and  relation  between  equivalent  and  effective 
descriptions 

It  follows  straightforwardly  from  (73)  that  the  dispersion  relation  p(&>)  may  be  re-written  as 


(f  =  ofpeqseq , 


(74) 


which,  after  using  Eqs.  (68)  and  (66),  is  proven  to  coincide  with  Eq.  (71). 

In  addition,  we  can  define  the  effective  characteristic  impedance  of  the  array  for  eigen-modal 
propagation: 


deff 


P!k0-c0 

(3lk0-c0 


(Zeff  +  Xeff) 


[Xeff  Xeff ) 


(75) 


In  the  absence  of  bianisotropic  effects  in  the  inclusion  aem  =  /f  =  0  and  the  last  side  in  Eq.  (75) 
simply  becomes 


deff 


(76) 


In  absence  of  bianisotropy,  therefore,  the  characteristic  impedance  is  not  directly  affected  by 
magneto  electric  coupling  at  the  lattice  level,  and  the  same  characteristic  impedance  is  obtained 
using  either  the  ratio  of  effective  or  equivalent  parameters.  Using  Eq.  (67),  we  may  write  in  the 
general  case 

Pav  _(%-goHfr+(4r+4r)  (?7) 

M av  iyMeff  ~  Pd  )  deff  ( Xeff  ~~  Xeff  ) 


which,  for  aem  =  /f  -  0 ,  may  be  also  written  as 


deff 


Peq  ~ 


(78) 


Eqs.  (74)  and  (78)  coincide  with  the  classic  retrieval  procedures  used  to  determine  effective 
permittivity  and  permeability  of  a  metamaterial  sample  from  their  secondary  parameters,  i.e., 
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their  eigen-modal  wave  number  (3  and  characteristic  impedance  rjr(;  [173],  [182],  This  implies 

that  the  equivalent  parameters  derived  in  (73)  from  first-principle  considerations  coincide  with 
those  retrieved  after  assuming  that  the  metamaterial  may  be  described  by  a  simple  model 
consisting  of  permittivity  and  permeability,  effectively  embedding  the  effects  of  %°eff  into  the 

equivalent  parameters.  In  souce-free  problems  in  which  the  sources  are  all  outside  the 
metamaterial  sample,  as  in  classic  retrieval  schemes,  it  is  tempting  to  put  aside  the  magneto¬ 
electric  coupling  coefficient  xljr ,  ar|d  use  the  equivalent  parameters  to  model  the  array 

scattering.  This  is  indeed  possible,  and  from  the  scattering  point  of  view  the  effective  and 
equivalent  description  is  perfectly  equivalent  in  this  source-free  scenario,  since  the  secondary 
parameters  of  the  supported  eigen-mode  coincide  in  the  two  descriptions.  However,  our  theory 
shows  that  the  equivalent  parameter  description  has  a  limited  physical  meaning  and  it  should  not 
be  used  to  characterize  the  electric  and  magnetic  response  of  the  metamaterial,  since  it  hides  an 
inherent  fonn  of  magneto-electric  coupling.  It  is  not  surprising  that  the  frequency  dispersion  of 
the  equivalent  parameters  may  contain  non-physical  artifacts  and  not  satisfy  passivity, 
reciprocity  or  causality  constraints  typical  of  local  parameters  [184], 

As  a  final  remark  with  respect  to  retrieval  procedures,  it  is  worth  noticing  that,  even  at 
frequencies  where  spatial  dispersion  and  magneto-electric  coupling  are  negligible  and  we  can 
write 


D  =s  E 

av  eq  a 


(79) 


as  in  a  natural  material,  our  theory  shows  that  the  averaged  fields  are  effectively  related  to  the 
microscopic  polarization  fields  through  Eq.  (56),  i.e.,  standard  retrieval  techniques  implicitly 
assume 


goE  +  Pj ■ 
E-Fh/s0 

AqH  + 

H-M  E/ Mo 


(80) 


In  other  words,  the  nature  of  the  averaged  polarization  currents  within  each  unit  cell,  whether 
electric  or  magnetic,  inherently  affects  the  definition  of  averaged  fields  used  to  calculate 
equivalent  or  effective  constitutive  parameters,  and  weak  spatial  dispersion  effects  associated 
with  artificial  magnetism  or  polarization  have  a  different  role  (contributing  to  Eav  and  Hav )  than 
the  direct  polarization  and  magnetization  vectors  (contributing  to  Dav  and  Bav ).  This  should  be 
considered  an  implicit  assumption  of  standard  retrieval  techniques  for  metamaterials. 


e.  Long-Wavelength  Limit  and  Relations  with  Approximate  Homogenization  Models 

For  sufficiently  long  wavelength  and  away  from  the  inclusion  resonances,  under  the  conditions 
k0d  <sc  1 ,  /3d  «;  1 ,  the  magneto-electric  coupling  effect  becomes  negligible,  implying  that 

°.^r=0.  (81) 
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Under  this  simple  condition,  and  in  the  absence  of  bianisotropic  effects  at  the  inclusion  level 
aem  =  0 ,  the  array  may  be  described  by  simple  constitutive  parameters 


£  rr  £  -  £ n 

eff  eq  0 


Meff  =  Meq  =  Mo 


1+ 
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-3  b 


ap  -c 
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int  J 

3  A’ 


1  +  -  _J 

V  ~  Cjnt  J 


(82) 


which  coincide  with  generalized  Clausius-Mossotti  relations  fully  taking  into  account  the 
coupling  among  the  inclusions  [176].  Under  condition  (81)  and  otem  =  0 ,  we  find  for  the 
eigenmodal  solution: 

p = 

P  £eff~£o  ’  (83> 

m  Meff -Mo 


which  is  exactly  coincident  with  Eqs.  (74)  and  (78)  for  equivalent  parameters,  ensuring  that  s 
and  jueq ,  as  well  as  the  retrieval  method  [173],  are  coincident  with  this  first-principle  definition 

of  effective  constitutive  parameters  in  the  absence  of  impressed  sources,  after  the  relevant 
assumption  that  magneto-electric  coupling  effects  stemming  from  spatial  dispersion  in  C'em  may 
be  neglected. 

If  also  Cint  may  be  considered  independent  of  [3  in  the  long-wavelength  limit,  then 


C„(ffl,/9->0)  =  -i5-  +  yi,  (84) 

3d  6k 

which  proves  that  Eqs.  (82)  and  this  homogenization  method  convege  to  local  classic  Clausius- 
Mossotti  formulas  for  periodic  arrays  [158]-[160]  when  co,f3  — >  0 .  In  the  following  numerical 
examples,  we  show  that  this  assumption  does  not  necessarily  hold  even  in  the  very  long 
wavelength  regime,  and  that  the  rigorous  homogenization  approach  introduced  here  may  provide 
results  significantly  different  from  quasi-static  approaches  even  for  kQd  <  1 ,  fid  <  1 . 


f.  Spatial  Dispersion  and  Extreme  Metamaterial  Parameters 

At  the  other  limit,  metamaterials  are  of  most  practical  interest  when  their  constitutive  parameters 
assume  extreme  (very  large,  very  low  or  negative)  values,  which  usually  arise  around  the 
inclusion  resonances.  The  homogenization  model  described  here  is  very  general,  and  in  principle 
applicable  to  any  value  of  co,fi .  However,  the  same  definition  of  homogenized  parameters,  as 
those  presented  in  the  previous  section,  inherently  neglects  the  array  granularity.  This  is 
particularly  relevant  near  these  resonances,  since,  despite  a  small  (k0d) ,  the  effective  eigen¬ 
modal  wavelength  may  become  comparable  with  the  period  as  (fid)  increases.  Although  these 
resonant  regions  are  quite  limited  in  bandwidth  for  passive  inclusions  in  the  long-wavelength 
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regime,  it  is  in  these  regions  that  the  effects  of  %°ff  and  of  inherent  spatial  dispersion  are  most 
relevant. 


i.  Near-zero  effective  material  parameters 

Limiting  ourselves  to  the  lossless  scenario  for  clarity,  consider  first  the  low-index  regions,  for 
which  (3d  -  0  for  finite  k0d  ,  of  interest  in  a  variety  of  applications  [  1 92]-[  1 97].  This  includes  s- 
near-zero,  u-near-zcro  and  low-index  metamaterials.  In  this  frequency  range,  (3  passes  from 
imaginary  to  real  values,  since  one  of  the  two  equivalent  parameters  crosses  the  real  axis.  For 
fi  =  0,  C'em  =  Xejj  -  0  and  their  value  also  goes  from  imaginary  to  real.  This  makes  sure  that  s 
and  p  stay  real  (and  one  of  them  negative)  in  Eq.  (73),  also  when  j3  e  3 .  As  shown  in  some  of 
the  following  numerical  examples,  crossing  this  zero-index  region  provides  significant 
deviations  between  the  equivalent  parameters  [seq,peq )  and  the  effective  parameters  [seff,peff} , 

as  a  symptom  of  inherent  spatial  dispersion,  consistent  with  the  results  in  [179].  It  should  be 
stressed  that,  surprisingly  enough,  in  this  region  j3,  C'em,xeff  are  all  very  close  to  zero,  implying 

very  long  effective  wavelengths  and  weak  magneto-electric  coupling;  however,  the  ratio  xlff  /  P 
(corresponding  to  Keff  in  the  notation  of  [184]),  crucial  for  determining  the  effect  of  spatial 

dispersion  in  the  denominator  of  (73),  is  not  necessarily  small,  providing  relevant  non-local 
effects  in  the  equivalent  parameters. 

ii.  Effective  parameters  near  the  bandgap  regions 

Another  region  of  interest  for  metamaterial  applications  is  the  one  near  the  edge  of  the  lattice 
bandgap,  for  which  (3d  —  n  .  Around  this  region,  large  positive  or  negative  values  of  permittivity 
and  penneability  may  be  obtained,  of  interest  in  a  variety  of  applications  [171],  [  1 98]-[  1 99] .  It  is 
evident  that  in  this  scenario  the  interaction  of  the  electromagnetic  wave  with  the  inclusions  may 
become  very  complex,  and  an  averaged  description  may  not  provide  much  insight  into  the 
physical  behavior  of  an  eigenmode  that  flips  its  phase  within  the  individual  unit  cell.  In 
particular,  inside  the  bandgap  the  same  definition  of  constitutive  parameters  is  not  meaningful, 
and  they  become  complex  even  for  lossless  inclusions.  It  is  relevant  to  study,  however,  the 
transition  between  the  homogenization  regime  and  the  bandgap  region,  where  extreme 
metamaterial  parameters  are  expected.  It  is  in  this  transition  region  that  this  rigorous 
homogenization  technique  becomes  particularly  important,  since  here  weak  spatial  dispersion 
effects  become  very  relevant,  even  in  the  long-wavelength  limit  k0d  <  1 .  Exactly  at  the  bandgap 

edge  the  periodic  properties  of  Cem  require  that  [176] 

Cem{P  =  7r  I  d)  =  0  Vco.  (85) 

For  center- symmetric  inclusions  ( aem  =0),  this  implies  that  the  general  dispersion  relation  (71) 
simplifies  into 

[or/1  -C(a>,  k ! -C(a>,  ;r/<7)]  =  0,  (86) 
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which  implies  that  a  bandgap  may  be  reached,  in  the  long-wavelength  limit  for  which  C  is 
small,  exclusively  near  an  electric  or  a  magnetic  resonance,  for  which  one  of  the  two  cT1  =  C  .  It 
follows  directly  from  (68),  (73)  that  at  such  resonance  one  of  the  equivalent  parameters 

Meg  =  Mo  (for  =  C)  or  seq  =  (for  aml  =  C  ),  (87) 


unless  the  two  resonances  are  degenerate,  as  considered  in  the  following  section. 
Correspondingly,  using  (74),  the  other  equivalent  parameter  has  to  become 


s 


eq 


(for  a;1 


C)  or  fueq  =/i0 


(for  a 


C )• 


(88) 


For  instance,  if  the  bandgap  associated  with  a  magnetic  resonance  is  considered,  as  it  is  the  case 
for  the  first  resonance  of  a  dielectric  inclusion  (see  example  2  in  the  following  numerical 
section),  the  eigen-wave  number  f5  and  the  corresponding  effective  penneability  jueff  rapidly 

increase  from  lower  frequencies  when  approaching  the  array  resonance,  until  they  hit  the 

n1 

bandgap  for  am  1  =  C  .  At  this  frequency  seq  =  s0 ,  /ueq  =  ju0 - j ,  and  therefore,  using  (73) 

(k0dy 


£eff  _  |  C0  Xeff 

£0  rcl[k0d) 

Meff  71  71  q 

_  /  .  /  /  J  C° 

ju0  k0a  (  k0a 


(89) 


Eq.  (75)  is  indeed  satisfied  by  Eq.  (89),  and  it  implies 

(90) 

770  k0  ykQd j 

which  suggests  that,  independent  on  the  geometry  of  the  inclusions,  at  the  bandgap  edge  the 
nonnalized  characteristic  impedance  simply  coincides  with  the  nonnalized  wave  number.  An 
analogous  derivation  for  bandgaps  at  an  electric  resonance  provides  the  inverse  of  Eq.  (90).  It  is 
evident  that  in  this  regime  /"ff  may  not  be  neglected  and  indeed  its  effect  is  comparable,  if  not 
larger,  than  the  one  of  seff  and  .  In  this  frequency  range  the  equivalent  parameters  (73)  lose 

their  physical  meaning  and  strongly  diverge  from  the  effective  parameters  (68),  as  discussed  in 
further  detail  in  [184]. 

It  is  evident  from  this  discussion  that  regions  with  extreme  (very  large,  very  low  or  negative) 
metamaterial  parameters  are  those  for  which  the  present  homogenization  technique  is  most 
relevant.  We  show  in  the  following  numerical  examples  how  this  rigorous  model  may  correctly 
model  the  exotic  features  of  metamaterials  within  these  frequency  bands,  capturing  the  weak 
spatial  dispersion  effects  that  are  usually  at  the  root  of  inconsistencies  in  approximate  and  less 
rigorous  homogenization  models. 


94 


Approved  for  public  release;  distribution  unlimited. 


g.  The  Special  Case  of  Matched  Inclusions  and  Degenerate  Resonances 

A  special  situation  of  interest  arises  in  the  case  in  which  the  electric  and  magnetic  resonances  of 
the  inclusions  are  degenerate,  as  it  may  happen  for  matched  inclusions,  which  support  combined 
electric -magnetic  resonances  at  the  same  frequency  [171], [173], [196].  These  are  particularly 
relevant  for  realizing  double  negative  metamaterials  with  matched  properties.  Under  this 
condition,  it  is  evident  that  ae  =  am  =  a  (at  least  over  a  given  range  of  frequencies  of  interest) 

and,  under  the  assumption  aem  =  0 ,  the  eigen-modal  dispersion  relation  (71)  simplifies  into 


Correspondingly,  Eq.  (72)  always  yields,  as  expectable 
m  !  P  =  do- 


The  effective  constitutive  parameters  are  simplified  into 
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p2  =kl*XL  =  kl^ 

£o  Ao 

deff  =  % 


(94) 


ensuring  that  effective  and  equivalent  parameters  are  matched,  as  expected. 

The  eigen-solution  f  of  (91)  is  necessarily  bounded  within  the  principal  period  |Re[/?]|  <n  /  d 

for  the  homogenization  to  hold.  In  the  lossless  limit  of  interest  here,  the  matched  condition 
ensures  that  there  are  no  finite  bandgaps  for  the  real  branches  of  ji ,  and  its  dispersion  bands  are 

connected  to  each  other  at  the  branchcuts  for  which  |Re[/?]|  =  n  /  d  .  The  absence  of  bandgaps  in 

the  eigen-modal  dispersion  is  particularly  interesting,  since  it  implies  that  the  effective 
parameters  are  always  real  defined,  producing  peculiar  dispersion  phenomena  in  this  special 
circumstance.  In  particular,  at  the  branch  cuts,  using  (85)  we  get 

al=C[co,j3-7r/d).  (95) 


Using  (93)  and  (66),  under  this  condition  we  find  the  interesting  relations 

8  eff  =  I1  eff  _  q 
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which  always  hold  at  branch-cut  points  |Re[/?]|  =  n  t  d  .  Notice  that  in  this  special  case  Eqs.  (87) 

-(88)  do  not  hold.  Although  the  effective  permittivity  and  penneability  are  zero,  the  magneto¬ 
electric  coefficient  is  large,  and  it  makes  sure  that  the  dispersion  relation  (94)  is  satisfied.  In  this 
regime,  by  locally  increasing  the  average  electric  field  would  curiously  affect  only  the  average 
magnetization,  and  dual  considerations  apply  to  the  average  magnetic  field.  It  is  evident  that  in 
this  scenario  the  effects  of  x°eff  anc*  of  weak  spatial  dispersion  cannot  be  neglected. 


Figure  44  -  (a)  Frequency  dispersion  of  the  electric  (thick  line)  and  magnetic  (thin)  normalized  polarizability  of  the 
individual  inclusions  for  the  four  metamaterial  arrays  considered  in  the  following  figures:  (solid)  dielectric  spheres 
with  permittivity  s  =  20  s0 ;  (dashed)  dielectric  spheres  with  s  =  120  sa ;  (dotted)  conducting  spheres;  (dash-dotted) 
magnetodielectric  spheres  with  s  =  20  s0  and  //  =  20  //„ ;  (b)  Ratio  of  electric  over  magnetic  polarizability  for  the 

same  geometries.  Flere  y  =  0.45  . 

Near  the  bandgaps,  there  are  two  finite  frequencies  satisfying  the  equation 
a'=C,„,±Cm,  (97) 

symmetrically  placed  around  the  branch-cut  point  aT1  =  C  (®,  n  /  tf) ,  at  which  all  the 

constitutive  parameters  diverge.  In  fact,  in  this  matched  scenario  the  effective  parameters  are  not 
bounded,  but 
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(98) 
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and  analogously  for  p  ,  ensuring  that  the  equivalent  parameters  in  (73)  remain  bounded  and 

satisfy  (94),  even  when  all  the  effective  constitutive  parameters  diverge.  These  two  closely 
spaced  vertical  asymptotes  for  the  effective  parameters  are  symmetrically  placed  around  any 
branch-cut  point  for  |Re[/?]|  =  n  /  d  and  they  are  separated  by  a  distinct  branch  that  crosses 

conditions  (96)  when  ax=C[co,P~nld\  Around  the  inclusion  resonance,  the  effective 
parameters  may  assume  zero  or  infinite  values  and  %°eff  cannot  be  neglected. 


Passed  any  of  these  resonances,  f3  assumes  negative  eigen-values,  which  corresponds  to  regions 
with  effective  negative  index  of  refraction,  in  which  both  the  effective  parameters  are  negative. 
At  the  connection  between  negative  and  positive  branches,  when  the  effective  parameters  cross 
zero,  in  this  matched  lossless  scenario  spatial  dispersion  may  be  safely  neglected,  since 
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Only  in  this  situation  we  may  achieve  a  truly  local  zero  pennittivity  and  penneability,  impedance 
matched  with  free-space,  useful  in  various  applications  envisioned  in  recent  papers  [196]-[197]. 
In  contrast,  in  the  unmatched  scenario,  as  outlined  above,  the  effects  of  spatial  dispersion  may 
not  be  neglected  in  regions  with  near-zero  effective  index  of  refraction. 


h.  Numerical  Examples  and  Further  Discussion 

In  this  section,  we  discuss  the  homogenization  of  four  specific  metamaterial  geometries  of 
interest.  Although  the  previous  formulation  is  applicable  to  arbitrary  lossy,  bianisotropic 
inclusions,  arbitrary  source  distribution  and  any  choice  of  ( co,  fi ) ,  here  we  focus  on 

metamaterials  composed  of  lossless  spheres  and  eigen-modal  propagation.  This  has  the 
advantage  of  providing  a  clearer  picture  of  the  difference  between  this  homogenization  approach 
and  other  quasi-static  techniques  and  retrieval  methods  focused  on  the  absence  of  embedded 
sources  in  the  array,  in  order  to  highlight  the  relevant  effects  of  spatial  dispersion  and  magneto¬ 
electric  coupling  coefficient  in  the  different  propagation  regimes  and  transition  regions  near  the 
metamaterial  resonances.  Our  numerical  examples  are  tailored  to  highlight  the  difference 
between  quasi-static  homogenization  models  and  this  rigorous  homogenization  approach.  For 
this  reason,  we  limit  our  analysis  to  dipolar  approximations  and  eigen-modal  propagation  in  the 
long-wavelength  regime,  of  major  interest  in  the  metamaterial  community.  In  future  works  we 
will  apply  the  general  multipolar  approach  introduced  in  Section  2  to  arbitrary  metamaterial 
inclusions  and  extend  our  numerical  analysis  to  the  presence  of  embedded  sources.  Since  we  deal 
with  spherical  particles,  we  can  use  analytical  closed-fonn  expressions  for  ae ,  am  [200].  Finally, 

for  the  reasons  outlined  above  and  the  purpose  of  this  paper,  we  focus  all  our  examples  in  the 
long- wavelength  region,  usually  considered  safe  for  quasi-static  homogenization  models,  i.e., 
(k0<7)<l  [162]-[164].  The  parameter  y  =  ald  is  also  introduced,  which  defines  the  ratio  of 
sphere  radius  over  lattice  period,  as  a  measure  of  the  array  density. 
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Figure  45  -  Frequency  dispersion  of  the  guided  wave  number,  and  its  approximations  as  defined  in  the  text,  for  an 
array  of  dielectric  spheres  with  s  =  20  s0 ,  with  (a)  y  =  0.45 ,  (b)  y  =  0.3  . 

Figure  44a  shows  in  logarithmic  scale  the  amplitude  of  the  normalized  polarizability  coefficients 
(thick  lines  for  the  normalized  electric  polarizability,  thin  for  the  magnetic  one)  for  four  different 
geometries  of  interest:  (1,  solid  lines)  dielectric  spheres  with  relative  pennittivity  sr  =20  and 
permeability  //,.  =  1 ;  (2,  dashed)  dielectric  spheres  with  sr  =  120  and  jur  =  1 ;  (3,  dotted)  perfectly 
conducting  spheres;  (4,  dash-dotted)  magneto-dielectric  spheres  with  sr  =  //.  =  20 .  For 
convenience,  Fig.  lb  shows  the  ratio  |«e|/|am|,  in  order  to  highlight  the  ratio  between  electric 
and  magnetic  response  at  the  inclusion  level.  Both  plots  show  their  variation  as  a  function  of 
[kQd) ,  for  the  density  factor  ;/  =  0.45  .  The  chosen  geometries  represent  specific  situations  of 

interest  in  common  metamaterial  arrays:  in  case  1  (solid  lines),  a  regular  array  of  dielectric 
spheres  is  considered,  far  from  their  individual  resonances,  but  still  with  a  good  contrast 
compared  to  the  background:  a  dominant  electric  polarization  is  expected  all  over  the  spectrum 
of  interest;  in  case  2  (dashed),  the  permittivity  is  increased  to  support  a  magnetic  and  an  electric 
resonance  within  the  frequency  band  of  interest,  in  analogy  with  established  designs  to  realize 
negative  metamaterial  parameters  [171]:  in  this  case,  more  interesting  features  are  expected  in 
the  metamaterial  response  near  these  resonances.  As  expected,  the  electric  response  is  dominant 
for  lower  frequencies,  but  the  first  resonance  is  magnetic.  In  case  3  (dotted),  conducting  particles 
are  considered,  for  which  electric  and  magnetic  responses  are  comparable,  and  for  lower 
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frequencies  the  electric  response  is  exactly  twice  the  magnetic  one;  in  case  4  (dash-dotted), 
impedance  matched  inclusions  are  employed,  supporting  two  combined  electric-magnetic 
resonances  within  the  band  of  interest.  It  is  noticed  that  in  all  these  examples,  lossless  conditions 
(70)  strictly  apply. 

Figure  45  a  shows  the  dispersion  of  normalized  eigen-wave  number  (effective  index  of  refraction) 
for  the  array  1  with  y  =  0.45  .  The  figure  compares  the  exact  eigen-modal  solution  of  the  guided 
wave  number  /?/  k0  (solid  line),  as  obtained  from  Eq.  (71),  with  various  approximate  solutions 
obtained  neglecting  spatial  dispersion  and  magneto-electric  effects,  as  follows:  the  dashed  line 
refers  to  the  dispersion  of  f3em ,  obtained  neglecting  the  magneto-electric  coupling  term  C'em  ,  as  in 

Eq.  (81);  the  dotted  line  shows  J3CM ,  which  in  addition  neglects  the  dispersion  effects  in  Cint , 
implying  C'em  =  0  and  Cmt  as  given  by  Eq.  (84),  coinciding  with  the  quasi-static  Clausius- 
Mossotti  homogenization  model;  the  dash-dotted  line  refers  to  /?, ,  obtained  neglecting  the 
magnetic  polarizability  effects  associated  with  the  magnetism  of  the  inclusions  (which  is  small  in 
this  geometry),  but  still  using  the  exact  Cint  expression;  finally  the  dash-dot-dot  line  refers  to 

Pe_CM ,  which  neglects  the  magnetic  effects  and  uses  Eq.  (84)  for  Cmt .  We  consider  all  these 
approximate  expressions  to  show  how  the  different  spatial  and  frequency  dispersion  terms, 
usually  neglected  in  quasi-static  homogenization  models,  affect  the  metamaterial 
homogenization  accuracy.  As  expected,  all  these  expressions  converge  to  the  same  quasi-static 
limit  when  (&>,/?)— >0,  but  the  approximate  expressions  start  deviating  from  the  exact 
expression  of  [3  for  relatively  low  values  of  k0d  .  In  particular,  it  is  surprising  to  notice  how,  by 
neglecting  the  magnetic  polarizability  of  the  particles,  which  in  this  case  is  orders  of  magnitude 
smaller  than  the  electric  one  (see  Fig.  lb),  the  dispersion  of  /?,  /  k0  diverges  quite  drastically 

from  the  exact  model,  implying  that  the  small  magnetism  in  these  dielectric  particles  should  not 
be  neglected,  as  one  may  be  tempted  to  do  looking  at  Fig.  lb.  The  effects  of  non-locality  and 
spatial  dispersion  in  Cint  start  playing  a  role  much  earlier  in  frequency  than  one  would  generally 

expect  for  such  simple  topology,  comparing  fiCM  with  J3  .  In  comparison,  the  magneto-electric 
coupling  effects  have  a  much  weaker  role,  and  start  being  relevant  only  around  k0d  -  1 .  Figure 
45b,  in  comparison,  shows  the  same  curves  calculated  for  a  less  dense  array,  with  y  =  0.3 .  As 
visible,  the  trend  of  the  various  curves  is  quite  similar,  although  the  effects  of  spatial  dispersion 
are  less  relevant  here,  as  the  interaction  among  inclusions  is  weaker.  In  particular,  the  magneto¬ 
electric  coupling  effects  associated  with  C'em  are  negligible  all  over  the  range  of  frequency  of 

interest,  since  f3em  practically  coincides  with  [3  in  this  less  dense  configuration.  Non-local 
effects  in  Cint  and  the  influence  of  the  small  magnetic  properties  of  the  inclusions  have  still  some 
relevant  effects  in  this  less  dense  scenario. 

Figure  46  shows  the  eigen-modal  dispersion  of  effective  constitutive  parameters  for  this  array  for 
y  =  0.45 .  The  top  panel  compares:  the  effective  permittivity  seff  (solid  black  line);  sem , 

calculated  after  neglecting  the  magneto-electric  coupling  coefficient  C'em ,  as  in  Eq.  (82) 
(dashed);  shc,  calculated  neglecting  also  the  effects  of  spatial  dispersion  in  Cinl ,  but  still 
considering  its  dependence  on  co  for  f3  =  0  (dash-dotted);  sCM ,  obtained  using  the  quasi-static 
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expression  for  Cint  given  in  (84)  (dotted),  which  coincides  with  the  Clausius-Mossotti  definition 
for  periodic  arrays  derived  in  [158];  finally  s  (solid  light  green),  defined  as  in  (73).  All  these 

expressions  yield  a  purely  real  permittivity,  as  expected  from  the  lossless  assumption.  However, 
sCM  rapidly  diverges  from  the  first-principle  value  of  permittivity  seff  .  The  value  of  seff  actually 

decreases  with  frequency  for  any  k0d  <  0.65 ,  due  to  the  effects  of  frequency  and  spatial 
dispersion  in  the  interaction  constants  for  dense  arrays.  By  neglecting  the  magneto-electric 
coupling,  as  in  sem  does  not  produce  any  sensible  difference  in  the  prediction  of  seff ,  but  the 

effects  of  spatial  and  frequency  dispersion  of  the  interaction  constants  are  quite  relevant,  as  seen 
analyzing  the  divergence  of  sloc  and  even  more  of  sCM  from  seff .  Finally,  the  divergence  of  s 

from  the  correct  value  of  seff  is  a  symptom  of  non-negligible  spatial  dispersion  and  magneto¬ 
electric  coupling  in  the  array,  which  are  evidently  not  negligible  in  such  dense  arrays. 

In  comparison,  the  permeability  is  accurately  predicted  by  all  the  approximate  models,  and  even 
the  local  or  Clausius-Mossotti  approximations  predict  extremely  well  its  weak  dispersion,  due  to 
the  significantly  lower  magnetic  response  of  the  spheres  all  over  the  frequency  range  of  interest. 
Interestingly,  only  n  shows  a  moderate  deviation  from  /ueff ,  which  highlights  how  the  effects 

of  x°eff  may  not  be  neglected  even  in  this  long-wavelength  regime.  Finally,  the  value  of  %°eff 

(bottom  panel)  becomes  relevant  only  towards  the  higher  end  of  this  frequency  range,  explaining 
the  divergence  of  effective  and  equivalent  parameters. 

Figure  47  calculates  the  secondary  effective  parameters  of  this  material,  obtained  using  the 
different  homogenization  models  of  Fig.  3.  In  particular,  Fig.  47a  compares  the  exact  value  of 
nonnalized  wave  number  J3  /  k0,  as  from  Fig.  45a,  with  the  approximate  values  /?  /  k0  =  > 

where  the  pedix  i  stands  for  any  of  the  approximate  models  used  in  Fig.  46.  This  plot  offers 
several  interesting  insights.  First  of  all,  it  is  noticed  that  /3ejJ  follows  extremely  well  the 

dispersion  of  (3em,  consistent  with  the  weak  effects  of  C'em  on  the  effective  permittivity. 
However,  both  curves  moderately  diverge  from  the  correct  value  fi  /  k0  in  the  range 
0.5  <k0d  <1,  confinning  that  the  effects  of  cannot  be  neglected  in  this  frequency  range. 
The  Clausius-Mossotti  model  J3CM  fails  even  more  substantially.  Fig.  47b  compares  the 
corresponding  values  of  effective  characteristic  impedance  iji  /  r/0  =  Jjui  /  sj  .  As  noticed  in  the 
previous  section,  %°eff  does  not  play  a  direct  role  in  the  impedance  when  xljj  =  0 ,  as  here,  and 
therefore  the  parameters  obtained  neglecting  C'em  yield  an  accurate  approximation  of  the 
effective  impedance  tjeff .  It  should  be  noted,  however,  that  the  relation  between  rjeff  and 

Pav  /  Mav  may  not  be  assumed  as  simple  as  (83),  due  to  the  effects  of  xljj  for  relatively  larger 

frequencies.  In  the  less  dense  array  case  of  Fig.  45b  (not  reported  here  for  brevity),  as  expectable 
the  effects  of  non-locality  and  spatial  dispersion  are  much  less  relevant,  but  still  Clausius- 
Mossotti  homogenization  formulas  would  deviate  considerably  from  the  proper  effective 
parameters.  The  use  of  equivalent  parameters,  that  embed  the  magneto-electric  effect  in  their 
same  definition  (73),  would  predict  correct  values  of  these  secondary  parameters,  consistent  with 
(74)-(75),  ensuring  that  their  use  for  scattering  purposes  in  the  eigen-modal  operation  is  perfectly 
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legitimate,  if  one  avoids  assigning  them  the  usual  physical  meaning  of  local  constitutive 
parameters. 


Figure  46  -  Frequency  dispersion  of  the  effective  constitutive  parameters,  and  their  approximations  as  defined  in  the 

text,  for  the  array  of  Fig.  45a. 
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Figure  47  -  Frequency  dispersion  of  the  effective  wave  number  and  characteristic  impedance  calculated  from  the 

constitutive  parameters  of  Fig.  46. 


Figure  48  -  Frequency  dispersion  of  the  guided  wave  number,  and  its  approximations  as  defined  in  the  text,  for  an 
array  of  dielectric  spheres  with  s  =  120  s0  and  y  =  0.45  .  The  thin  solid  line  corresponds  to  the  imaginary  branch  of 

P. 
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Figure  49  -  Frequency  dispersion  of  the  effective  constitutive  parameters,  and  their  approximations  as  defined  in  the 
text,  for  the  array  of  Fig.  5.  Dashed  lines  in  the  bottom  panel  refer  to  branches  with  imaginary  values. 

Consider  now  the  second  metamaterial  of  interest,  composed  of  spheres  with  sr  =120,  which 
support  a  magnetic  and  an  electric  resonance  within  the  low  frequency  regime  considered  here. 
Figure  48  shows  the  eigen-wave  number  dispersion  for  such  array  with  y  =  0.45  ,  with  symbols 
analogous  to  Fig.  2.  It  is  immediately  recognized  that  the  exact  dispersion  of  normalized  wave 
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number  [3  /  k0  (solid  lines)  is  much  more  intricate  than  in  the  previous  situation.  As  expected, 
/3  /  k0  initially  grows  with  frequency,  until  hitting  the  first  band-gap  of  the  array  at 
( k0d )  =  0.594  ,  at  the  magnetic  resonance  frequency  aj  -C  =  0  .  The  narrow  frequency  region 

within  the  bandgap  should  be  completely  disregarded  in  terms  of  homogenization,  since,  as 
discussed  above,  the  effects  of  array  granularity  plays  a  major  role,  and  the  same  definition  of 
averaged  constitutive  parameters  loses  its  macroscopic  meaning.  Passed  the  magnetic  bandgap,  a 
branch  with  imaginary  wave  number  /3  =  -  //?  is  reached  (thin  solid  line),  which  connects  with 

the  next  real  branch  at  ( k0d )  =  0.723 ,  at  the  point  for  which  J3  =  0  .  The  following  bandgap  is 
then  hit  at  the  electric  resonance  frequency  a~p  -  C  =  0 ,  at  (kltd )  =  0.891 ,  and  the  next  real 
branch  is  obtained  at  (k0d)  -  0.909 .  As  seen  in  Fig.  7,  this  behavior  is  well  described  by 
approximate  dispersion  relations,  even  neglecting  the  effects  of  C'em  or  even  the  spatial 
dispersion  in  Cint ,  as  in  f3em  and  /3CU  respectively,  since  the  local  inclusion  resonances  now 
dominate  the  spatial  dispersion  effects  of  the  lattice.  Of  course,  in  this  scenario  it  is  not  possible 
to  neglect  the  magnetic  effects  in  the  dielectric  particles,  as  for  (3e  and  f3e  CM ,  since  this  would 
completely  miss  the  first  magnetic  bandgap  resonance. 

The  effective  constitutive  parameters  of  this  array  are  reported  in  Fig.  49,  with  analogous 
symbols  as  in  the  previous  example.  If  the  spatial  dispersion  effects  are  negligible  in  evaluating 
p(co)  in  Fig.  48,  they  play  a  major  role  in  the  correct  definition  of  the  constitutive  parameters,  in 

particular  near  the  electric  and  magnetic  resonances  of  the  inclusions.  First,  it  is  noticed  that 
Clausius-Mossotti  fonnulas  completely  miss  the  relevant  magneto-electric  coupling  arising  near 
the  bandgaps,  and  the  pennittivity  especially  suffers  of  this  approximation,  starting  from  very 
low  frequencies.  Towards  the  first  (magnetic)  resonance,  sem  may  approximate  relatively  well 

the  effective  pennittivity  seff,  confirming  that  the  effect  of  C'em  is  small  on  the  permittivity 

dispersion,  dominated  by  the  local  inclusion  resonances.  However,  the  value  of  xlff  assumes 

large  values  near  the  two  resonances  and  it  cannot  be  neglected  in  a  proper  homogenization 
model.  Near  the  magnetic  resonance,  the  effective  pennittivity  experiences  a  sharp  Lorentzian 
resonance,  completely  missed  by  sCM  and  even  by  sloc,  which  is  an  evident  symptom  of  the 
magneto-electric  coupling  in  the  anay.  In  contrast,  the  various  models  for  magnetic  permeability 
all  have  good  agreements  with  the  effective  model  (with  the  exception  of  a  small  resonant 
feature  arising  at  the  electric  bandgap  resonance  of  the  array).  In  the  region  where  (3  is 
imaginary,  immediately  following  the  bandgaps,  all  the  models  correctly  predict  a  negative 
effective  penneability  or  permittivity  region,  which  crosses  zero  at  (k0d)  -  0.723  and 

( k0d )  =  0.909  ,  together  with  the  value  of  /A  In  this  negative  parameter  range,  as  expected, 

is  imaginary  (dashed  lines  in  the  bottom  panel),  which  ensures  that  the  equivalent  parameters  are 
real  quantities  (one  of  them  negative). 

Special  attention  should  be  paid  to  the  dispersion  of  the  equivalent  permittivity  s  in  Fig.  51a 

(lighter  green  line).  Its  slope  is  negative  all  the  way  until  the  magnetic  bandgap,  producing  an 
anti-resonance  effect  similar  to  the  one  extracted  from  retrieval  procedures  in  various 
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metamaterial  geometries  near  magnetic  resonances  [173]-[1 82] .  It  is  evident  from  this  analysis 
that  these  effects  are  associated  to  incorporating  x'\j  in  the  definition  of  equivalent  pennittivity. 


Figure  50  -  (Color  online):  Frequency  dispersion  of  the  effective  wave  number  and  characteristic  impedance 

calculated  from  the  constitutive  parameters  of  Fig.  49. 

It  is  true  that  the  two  equivalent  parameters  may  describe  well  the  guidance  properties  of  the 
array,  but  their  physical  meaning  in  this  case  diverges  considerably  from  the  first-principle 
definition  of  permittivity  and  penneability.  A  simple  metamaterial  model  based  on  just 
equivalent  permittivity  and  permeability  would  fail  capturing  the  physics  of  the  array  near  the 
first  bandgap  resonance,  predicting  s  =  £0  (87),  when  in  reality  the  averaged  polarization 

vector  has  a  strong  resonance  cancelled  by  the  coupled  magnetization.  These  effects  have  been 
discussed  in  more  detail  in  [184].  If  the  discrepancy  between  s  and  seff  may  have  been 

expected  near  this  resonant  regime,  another  transition  region  in  which  the  equivalent  parameters 
seq ,  /ueq  lose  their  usual  meaning  is  the  region  near  (  k()d)  -  0.723 ,  for  which  [3  —  0 .  As 
confirmed  by  Fig.  8a,  and  anticipated  above,  in  this  region: 
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(100) 


S  it  —  S  —  Si 

ejj  em  loc 

Heff  ~  0 

Indeed,  the  correct  value  of  effective  permittivity  coincides  with  the  local  value  sloc,  since 
P  -  0  ,  but  this  value  is  substantially  different  from  s  .  This  is  due  to  the  fact  that,  although  the 
magneto-electric  coefficient  is  near  zero,  the  ratio  %°ff  /  p  is  finite,  implying  that  seq  diverges 
from  eeff  and  it  loses  the  meaning  of  the  average  electric  polarizability  of  the  array.  In  this  near¬ 
zero  index  region,  the  weak  spatial  dispersion  represented  by  %°eff  may  not  be  neglected, 

although  the  effective  wavelength  is  very  large.  This  confirms  the  results  in  [179]  derived  for 
periodic  arrays  of  split-ring  resonators,  which  discussed  the  presence  of  non-negligible  spatial 
dispersion  effects  in  this  extremely  long-wavelength  (pd  «;  1 )  scenario. 


Figure  51  -  Frequency  dispersion  of  the  guided  wave  number,  and  its  approximations  as  defined  in  the  text,  with 
frequency  for  an  array  of  conducting  spheres,  for  (a)  y  =  0.45  ,  (b)  y  =  0.3  . 
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Figure  52  -  Frequency  dispersion  of  the  effective  constitutive  parameters,  and  their  approximations  as  defined  in  the 

text,  for  the  array  of  Fig.  8a. 
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Figure  53  -  Frequency  dispersion  of  the  effective  wave  number  and  characteristic  impedance  calculated  from  the 

constitutive  parameters  of  Fig.  52. 

Figure  50  shows  the  dispersion  of  the  effective  index  of  refraction  and  characteristic  impedance 
obtained  through  the  various  parameters  in  Fig.  49,  analogous  to  Fig.  47.  As  it  is  visible,  all  the 
curves  agree  with  high  accuracy  within  the  real  branches,  since  their  dispersion  is  dominated  by 
the  local  resonances  at  the  inclusion  level.  This  example  clearly  shows  that  indeed  f5  and  ij  for 
this  array  may  be  easily  derived  applying  local  concepts,  like  Clausius-Mossotti  relations  or 
simple  retrieval  procedures,  since  they  are  dominated  by  local  resonances  at  the  inclusion  level; 
however,  inferring  from  these  secondary  parameters  the  physical  values  of  pennittivity  and 
penneability,  as  commonly  done,  may  lead  to  physical  artifacts  and  inconsistencies  [184],  since 
in  the  resonant  regions  the  effects  of  x"u  cannot  be  neglected. 

As  a  third  example,  consider  the  case  of  an  array  of  conducting  particles,  as  in  Fig.  44,  dotted 
lines.  Figure  51  shows  the  dispersion  of  wave  numbers  for  y  =  0.45  (a)  and  7  =  0.3  (b), 
analogous  to  Fig.  45.  In  this  case,  the  wave  numbers  predicted  using  just  electric  effects  of  the 
particles  are  evidently  incorrect,  since  the  magnetic  contribution  for  conducting  particles  is  never 
negligible.  Moreover,  the  effect  of  the  coupling  coefficient  C'em  is  particularly  relevant  in  this 
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conducting  scenario,  which  shows  the  maximum  divergence  between  (3  and  (3 em ,  due  to  the 
relevance  of  the  magnetic  effects  even  at  very  low  frequencies. 

Figure  52  shows  the  corresponding  constitutive  parameters  for  the  case  y  =  0.45  .  seff  also  in  this 

scenario  shows  a  distinct  negative  slope,  all  over  the  range  of  frequencies  considered  here,  due  to 
non-negligible  spatial  and  frequency  dispersion.  This  is  compensated  by  the  positive  slope  of  the 
effective  permeability,  which  assumes,  as  expected,  diamagnetic  values  [201],  Only  the 
Clausius-Mossotti  quasi-static  model  predicts  a  positive  dispersion  for  sCM  ,  whereas  all  the  other 

models  correctly  follow  the  trend  of  seff .  As  good  confirmation  of  the  strong  influence  of  xljf , 
the  deviation  of  the  equivalent  parameters  s  ,  and  in  particular  jueq ,  from  the  effective 

parameters  is  pretty  relevant.  Figure  10,  finally,  shows  the  dispersion  of  the  calculated  wave 
numbers  and  characteristic  impedances  obtained  using  the  effective  constitutive  parameters  of 
Fig.  52.  It  is  seen  how  all  the  curves  agree  reasonably  well  with  the  exact  dispersion  of  rje(f , 

except  the  quasi-static  Clausius-Mossotti  formula,  which  neglects  frequency  and  spatial 
dispersion  effects  of  the  interaction  constants.  The  divergence  of  all  the  curves  from  the  exact 
dispersion  of  (3  is  particularly  striking,  as  a  symptom  of  the  relevance  of  the  magneto-electric 

coefficient  xlg-  •  We  have  also  analyzed  the  less  dense  configuration  with  ;/  =  0.3  ,  as  in  Fig.  51b 

(not  reported  here  for  brevity),  which  indeed  provides  analogous  results,  but  less  strong 
variations  from  the  background  parameters,  as  expected. 


K)d 


Figure  54  -  Frequency  dispersion  of  the  guided  wave  number,  and  its  approximations  as  defined  in  the  text,  for  an 
array  of  magnetodielectric  spheres  with  s  =  20 s0 ,  fu  =  2Qju0,  y  =  0.45  . 


Finally,  as  the  fourth  relevant  metamaterial  geometry  of  interest,  Figs.  54-55  show  the  case  of  an 
array  of  magnetodielectric  matched  spheres,  as  in  Fig.  1  (blue  lines),  which  are  able  to  produce 
degenerate  electric-magnetic  resonances  at  the  same  frequency,  leading  to  negative  index  bands 
of  operation  [171],  [173].  Figure  54  shows  the  eigen- wave  number  dispersion  for  this  array, 
consistent  with  the  previous  examples.  Indeed,  the  wave  number  experiences  resonances  near  the 
combined  electric-magnetic  resonances  of  the  inclusions,  as  predicted  by  (91).  In  particular,  /? 
grows  from  its  quasi-static  value  until  reaching  the  first  branch-cut  at  f}  =  n  !  d  ,  at  the  resonant 
condition  (95).  Exactly  as  formulated  in  the  previous  section,  the  effective  constitutive 
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parameters  (Fig.  55)  diverge  twice  around  this  resonance  condition,  at  the  frequencies  (91),  and 
the  corresponding  wave  number  enters  a  region  of  negative  effective  index  of  refraction.  As 
expected,  also  x'ejj  diverges  near  this  resonance,  and  it  cannot  be  neglected  in  a  proper 

homogenization  model.  All  the  effective  parameters  are  real  for  any  frequency,  and  their  value  is 
unbounded,  as  described  in  the  previous  section.  Consistently,  the  approximate  models  of  [i  (in 
this  scenario  they  are  identically  equal  to  the  dispersion  of  the  permittivity  models  in  Fig.  55a, 
since  permittivity  and  permeability  are  matched)  may  considerably  fail  near  these  resonances.  As 
already  mentioned  in  the  previous  section,  in  this  matched  geometry  the  eigen-modal 
characteristic  impedance  is  always  r/0 ,  equal  to  the  ratio  Mav  /  Pav .  Therefore,  classic  retrieval 
techniques,  as  well  as  the  one  employed  in  [173],  may  derive  with  good  approximation  the 
behavior  of  seff,  jue/f  far  from  the  branch-cuts,  for  which  and  the  effects  of  spatial 

dispersion  are  minimal.  Flowever,  near  the  branch-cut  resonances  a  more  refined  model  as  the 
one  described  here  is  required  for  a  proper  physical  description  of  the  array. 


Figure  55  -Frequency  dispersion  of:  (a)  effective  permittivity  and  its  approximations  (the  normalized  permeability 
and  effective  wave  number  have  exactly  the  same  dispersion  in  this  case)  and  (b)  magneto-electric  coefficient  for  the 

array  of  Fig.  54. 
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i.  Conclusions 


We  have  laid  out  in  this  section  the  foundations  of  a  first-principle  homogenization  theory  to 
properly  define  quasi-local  effective  constitutive  parameters  for  a  periodic  array  of 
magnetodielectric  and  conducting  inclusions.  We  have  shown  that  a  proper  and  rigorous 
description  of  the  wave  interaction  in  metamaterial  arrays  requires  taking  into  account  weak 
spatial  dispersion  effects  and  magneto-electric  coupling  within  the  lattice,  even  in  the  long- 
wavelength  limit  and  for  center-symmetric  inclusions.  In  this  context,  we  have  derived  closed- 
form  analytical  formulas  for  the  rigorously  defined  effective  constitutive  relations  of  an  arbitrary 
metamaterial  array,  highlighting  the  limitations  of  commonly  used  models  that  neglect  these 
effects.  Although  our  theory  is  very  general,  the  numerical  results  presented  here  are  focused  on 
isotropic  metamaterial  arrays,  lossless  materials  and  eigen-modal  propagation,  in  order  to  better 
highlight  the  specific  effects  of  spatial  dispersion,  above  all  in  the  transition  between  different 
eigen-modal  regions,  in  the  long-wavelength  limit.  We  have  highlighted  in  our  calculations  the 
effects  of  neglecting  spatial  and  frequency  dispersion  in  the  interaction  constants  and  in  the 
constitutive  models  of  metamaterials  formed  by  magnetodielectric  and  conducting  spheres, 
which  may  cause  severe  artifacts  and  divergence  from  the  first-principle  effective  parameters  of 
the  metamaterial.  We  believe  that  these  results  may  be  fundamental  for  a  proper  description  of 
artificial  materials  with  nonconventional  properties  in  complex  environments,  for  their  realistic 
applications  in  practical  devices  of  interest  to  the  U.S.  Air  Force. 


8.  General  Conclusions 

We  have  reported  here  a  series  of  basic  research  results  as  the  outcome  of  our  extensive  efforts 
sponsored  by  the  U.S.  Air  Force  Research  Laboratory  with  Contract  No.  FA8718-09-C-0061. 
We  have  focused  in  this  report  on  the  most  relevant  aspects  of  our  findings,  and  in  particular  on 
guided  and  leaky-wave  modes  along  linear  arrays  of  nanoparticles,  also  considering  and 
modeling  the  realistic  presence  of  technological  disorder,  comparison  of  the  guidance  properties 
along  linear  and  planar  arrays  of  nanoparticles  and  nanovoids  in  different  realistic  geometries, 
guided  and  leaky  modes  along  parallel  arrays  of  nanoparticles,  propagation  along  periodic  arrays 
of  nanoparticles  and  their  rigorous  homogenization  as  metamaterials  and  nonconventional 
artificial  materials.  We  believe  that  these  findings  may  be  of  significant  importance  to  realize 
novel  optical  devices  and  exotic  materials  for  a  variety  of  applications  of  interest  to  the  U.S.  Air 
Force.  A  complete  overview  of  our  results  may  be  found  in  the  series  of  papers  we  have 
published  during  this  effort,  as  reported  in  the  following  section. 
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